## Replication code for
## Linking Perspectives: A Field Experiment on the Role of Multi-Layer
## Networks in Refugee Information Sharing
## Clark, Larson and Lewis
## Created: Aug 28, 2023
## Updated: Apr 19, 2024



library(igraph)
library(xtable)
library(tidyverse)
library(stargazer)
library(broom)


## Load the replication data
D <- read_csv("ClarkLarsonLewis2024.csv") #1604 x 54

#########################################################################
## The next section of code creates household networks, one for each village


######
## Create village networks

## First, need to make sure the deidentified name variables which are entered
## as random numbers that stand for names are understood to be characters 
## and not numeric

## Household names:

D$ego <- as.character(D$ego)
D$hhh <- as.character(D$hhh)
D$hh1 <- as.character(D$hh1)
D$hh2 <- as.character(D$hh2)
D$hh3 <- as.character(D$hh3)
D$hh4 <- as.character(D$hh4)
D$hh5 <- as.character(D$hh5)
D$hh6 <- as.character(D$hh6)
D$hh7 <- as.character(D$hh7)
D$hh8 <- as.character(D$hh8)
D$hh9 <- as.character(D$hh9)
D$hh10 <- as.character(D$hh10)
D$hh11 <- as.character(D$hh11)

## Visit layer names:

D$visit1 <- as.character(D$visit1)
D$visit2 <- as.character(D$visit2)
D$visit3 <- as.character(D$visit3)
D$visit4 <- as.character(D$visit4)
D$visit5 <- as.character(D$visit5)

## Meal layer names:

D$meal1 <- as.character(D$meal1)
D$meal2 <- as.character(D$meal2)
D$meal3 <- as.character(D$meal3)
D$meal4 <- as.character(D$meal4)
D$meal5 <- as.character(D$meal5)

## Borrow layer names:

D$borrow1 <- as.character(D$borrow1)
D$borrow2 <- as.character(D$borrow2)
D$borrow3 <- as.character(D$borrow3)
D$borrow4 <- as.character(D$borrow4)
D$borrow5 <- as.character(D$borrow5)

## Rumor layer names:

D$rumor1 <- as.character(D$rumor1)
D$rumor2 <- as.character(D$rumor2)
D$rumor3 <- as.character(D$rumor3)
D$rumor4 <- as.character(D$rumor4)
D$rumor5 <- as.character(D$rumor5)

## And the names of refugee discussion partners offered

D$spokeref1 <- as.character(D$spokeref1)
D$spokeref2 <- as.character(D$spokeref2)
D$spokeref3 <- as.character(D$spokeref3)
D$spokeref4 <- as.character(D$spokeref4)
D$spokeref5 <- as.character(D$spokeref5)
D$spokeref6 <- as.character(D$spokeref6)
D$spokeref7 <- as.character(D$spokeref7)
D$spokeref8 <- as.character(D$spokeref8)
D$spokeref9 <- as.character(D$spokeref9)
D$spokeref10 <- as.character(D$spokeref10)


## Next, create separate village datasets since we'll
## build separate networks for each village

v1 <- D[which(D$village_number == 1),]
v2 <- D[which(D$village_number == 2),]
v3 <- D[which(D$village_number == 3),]
v4 <- D[which(D$village_number == 4),]
v5 <- D[which(D$village_number == 5),]
v6 <- D[which(D$village_number == 6),]
v7 <- D[which(D$village_number == 7),]
v8 <- D[which(D$village_number == 8),]
v9 <- D[which(D$village_number == 9),]
v10 <- D[which(D$village_number == 10),]
v11 <- D[which(D$village_number == 11),]
v12 <- D[which(D$village_number == 12),]


## Next, let's write a function which will take a vector of (deidentified) names
## that will serve as our egos and a list of other names that will be the alters
## and will return the two-column matrix that will be the edgelist for this network

make_el <- function(data, ego, alter_list){ #ego and alter_list both include the dataset too
  
  el <- matrix(data = NA, nrow = nrow(data)*length(alter_list), ncol = 2) #make empty 2-column matrix
  colnames(el) <- c("ego", "alter") #name the columns
  counter <- 1 #start a counter at 1
  
  for(i in 1: nrow(data)){
    for(j in 1:length(alter_list)){ #for all the alters listed in the list
      el[counter,] <- c(ego[i], alter_list[[j]][i]) #add a link between the ego and that alter
      counter <- counter + 1 #and advance the counter by 1 so we know where to put the next one
    }
  }
  
  notempty <- which(el[,2] != "") ## Now clean it so alters of "" are listwise deleted.
  el <- el[notempty,]
  
  return(el) 
}

## Next, let's write a function that takes an individual-level edgelist
## and converts it to a household-level edgelist.
## Works by taking a list of everyone in each household, comparing
## the individuals in the individual-level edgelist to it, and makeing a new
## edgelist with all links between households that the individual edgelist
# and household list imply.  So a link will be added between two households
## for every instance that someone in one household listed someone the other household.

convert_to_hh <- function(hhlist, el){ #takes a hh reference list and an indvidual edgelist to make hh nw
  
  hhel <- matrix(data = NA, nrow = 0, ncol = 2)
  
  #For the ego and then alter, identify the household(s) they're in
  ## in terms of the hh number in the hh list
  ## (this assigns node numbers in order that hh appears in the data; preserves data ordering)
  for(i in 1:nrow(el)){
    
    hhego <- c()
    for(j in 1:length(hhlist)){
      if(sum(el[i,1] == hhlist[[j]])>0){
        hhego <- c(hhego, j)
      }
    }
    
    hhalter <- c()
    for(k in 1:length(hhlist)){
      if(sum(el[i,2] == hhlist[[k]])>0){
        hhalter <- c(hhalter, k)
      }
    }
    
    # If the ego or alter wasn't in any of the households in hhlist, go to the next row in el
    ## (Automatically closes this network)
    if(is.null(hhego)){
      next
    }
    if(is.null(hhalter)){
      next
    }
    
    #Now add all household links implied by all pairs of ego + alter to hhel
    
    for(l in 1:length(hhego)){
      for(m in 1:length(hhalter)){
        hhel <- rbind(hhel, c(hhego[l], hhalter[m]))
      }
    }
    
    #print(i)  
  }
  
  return(hhel)
  
}

## Finally, let's write a function that takes a network and a dataframe
## and adds some certain node attributes to the network.  

addatts <- function(nw, data){
  
  V(nw)$trt <- data$trt
  V(nw)$wasrefugee <- data$been_refugee
  V(nw)$metrefugee <- data$refugee_meet
  V(nw)$threat1 <- data$threat1
 # V(nw)$religion <- data$relig
  V(nw)$language <-data$resp_language
  V(nw)$salient <- ifelse(data$important_surr_refug== 1, 1, 0)
  
  out <- nw
  
  return(out)
  
}

## Now we'll put all three functions together in a grand, network-builder function.
## This function calls the earlier functions to start from a dataframe
## and build all the networks we have data on at once: visit, meal, borrow, rumor,
## the union of these four, the network of whom they said they
## planned to tell (willtell), and finally the network of who spoke about refugees
## with whom in the endline (spokeref, sometimes we call it didtell).

## This function is flexible; it offers the options to make the networks directed or
## undirected network, to remove links that are duplicated, and to remove loops 
## (when someone lists someone in their own household).  These default to making
## an undirected network without duplicated links or loops

## Recall that this function builds the network in the order the nodes
## appear in the data.
## NOTE: THIS FUNCTION DIFFERS SLIGHTLY FROM THE VERSION IN OTHER CODE FILES.
## SPECIFICALLY, IT ASSIGNS THE EGO'S NAME TO AN ATTRIBUTE CALLED EGO
## AND ASSIGNS THE HOUSEHOLD NUMBER FROM THE DATA TO THE NODE'S NAME
makehhnws <- function(data, directed = FALSE, rmduplinks = TRUE, rmloops = TRUE){
  
  ## First grab the relevant variables for each network from the dataset
  visitlist <- list(data$visit1, data$visit2, data$visit3, data$visit4, data$visit5) 
  meallist <- list(data$meal1, data$meal2, data$meal3, data$meal4, data$meal5)
  borrowlist <- list(data$borrow1, data$borrow2, data$borrow3, data$borrow4, data$borrow5)
  rumorlist <- list(data$rumor1, data$rumor2, data$rumor3, data$rumor4, data$rumor5)
  unionlist <- list(data$visit1, data$visit2, data$visit3, data$visit4, data$visit5, 
                    data$meal1, data$meal2, data$meal3, data$meal4, data$meal5, 
                    data$borrow1, data$borrow2, data$borrow3, data$borrow4, data$borrow5,
                    data$rumor1, data$rumor2, data$rumor3, data$rumor4, data$rumor5)

  spokereflist <- list(data$spokeref1, data$spokeref2, data$spokeref3, data$spokeref4, data$spokeref5,
                       data$spokeref6, data$spokeref7, data$spokeref8, data$spokeref9, data$spokeref10)
  
  
  ## Now make the edgelists, individual, using the function built above called make_el
  visit <- make_el(data = data, ego = data$ego, alter_list = visitlist) #326
  meal <- make_el(data = data, ego = data$ego, alter_list = meallist) #219
  borrow <- make_el(data = data, ego = data$ego, alter_list = borrowlist) #147
  rumor <- make_el(data = data, ego = data$ego, alter_list = rumorlist) #143
  union <- make_el(data = data, ego = data$ego, alter_list = unionlist) #4437

  spokeref <- make_el(data = data, ego = data$ego, alter_list = spokereflist)#768
  
  ## Now convert the individual edgelists to household edgelists.  First step is
  ## to create a matrix of household members, and from that matrix, a list.
  
  ## First make a matrix of household members,
  familymat <- cbind(data$ego, data$hhh, data$hh1, data$hh2, data$hh3, data$hh4, data$hh5,
                     data$hh6, data$hh7, data$hh8, data$hh9, data$hh10, data$hh11)
  
  ## Now construct a list of household members from this matrix
  
  hhlist <- list()
  hhcount <- c()
  
  for(i in 1:length(data$ego)){
    #Make sure to match to names, not blank strings
    names <- which(familymat[i,] != "")
    #Find the names of family members for 
    membs <- familymat[i,names]
    
    hhlist[[i]] <- unique(membs) #grab the unique names because ego and hhh might be repeated
    hhcount[i] <- length(unique(membs))
  }
  
  ## Now convert each edgelist to household edgelist using the function built
  ## above called convert_to_hh by comparing the individual edgelist to the 
  ## household member list
  visithh <- convert_to_hh(hhlist, visit)
  mealhh <- convert_to_hh(hhlist, meal)
  borrowhh <- convert_to_hh(hhlist, borrow)
  rumorhh <- convert_to_hh(hhlist, rumor)
  unionhh <- convert_to_hh(hhlist, union)
  spokerefhh <- convert_to_hh(hhlist, spokeref)
  
  if(rmduplinks == TRUE){ #Purge duplicated links
    
    if(directed == TRUE){
      visithh <- unique(visithh, MARGIN = 1)
      mealhh <- unique(mealhh, MARGIN = 1)
      borrowhh <- unique(borrowhh, MARGIN = 1)
      rumorhh <- unique(rumorhh, MARGIN = 1)
      unionhh <- unique(unionhh, MARGIN = 1)
      spokerefhh <- unique(spokerefhh, MARGIN = 1)
    }
    
    ## Removing duplicate links from undirected network requires first ordering
    ## the egos and alters and then taking the unique
    ## Sorting each row of nxm returns mxn, so take transpose to return to nxm
    
    if(directed == FALSE){
      visithh <- t(apply(visithh, MARGIN = 1, sort))
      visithh <- unique(visithh, MARGIN = 1)
      mealhh <- t(apply(mealhh, MARGIN = 1, sort))
      mealhh <- unique(mealhh, MARGIN = 1)
      borrowhh <- t(apply(borrowhh, MARGIN = 1, sort))
      borrowhh <- unique(borrowhh, MARGIN = 1)
      rumorhh <- t(apply(rumorhh, MARGIN = 1, sort))
      rumorhh <- unique(rumorhh, MARGIN = 1)
      unionhh <- t(apply(unionhh, MARGIN = 1, sort))
      unionhh <- unique(unionhh, MARGIN = 1)
      spokerefhh <- t(apply(spokerefhh, MARGIN = 1, sort))
      spokerefhh <- unique(spokerefhh, MARGIN = 1)
    }
    
  }
  
  ## also deal with and count the loops
  
  if(rmloops == TRUE){ #purge loops
    vloops <- which(visithh[,1] == visithh[,2])
    if(length(vloops)>0){
      visithh <- visithh[-(vloops),]
    }
    
    mloops <- which(mealhh[,1] == mealhh[,2])
    if(length(mloops)>0){
      mealhh <- mealhh[-(mloops),]
    }
    
    bloops <- which(borrowhh[,1] == borrowhh[,2])
    if(length(bloops)>0){
      borrowhh <- borrowhh[-(bloops),]
    }
    
    rloops <- which(rumorhh[,1] == rumorhh[,2])
    if(length(rloops)>0){
      rumorhh <- rumorhh[-(rloops),]
    }
    
    uloops <- which(unionhh[,1] == unionhh[,2])
    if(length(uloops)>0){
      unionhh <- unionhh[-(uloops),]
    }
    
    sloops <- which(spokerefhh[,1] == spokerefhh[,2])
    if(length(sloops)>0){
      spokerefhh <- spokerefhh[-(sloops),]
    }
    
  }
  
  ## Now build the networks using the igraph function graph_from_edgelist
  visitnw <- graph_from_edgelist(visithh, directed = directed)
  mealnw <- graph_from_edgelist(mealhh, directed = directed)
  borrownw <- graph_from_edgelist(borrowhh, directed = directed)
  rumornw <- graph_from_edgelist(rumorhh, directed = directed)
  unionnw <- graph_from_edgelist(unionhh, directed = directed)
  spokerefnw <- graph_from_edgelist(spokerefhh, directed = directed)
  
  
  ## Add singletons that have a node index higher than the last included
  ## as a party to a link in the el
  ## (igraph fills in intermediate singletons, but not those that happen to
  ## be at the end of the node list)
  
  if(vcount(visitnw) < length(data$ego)){
    lastone <- max(V(visitnw)) #find the last vertex added
    numtoadd <- length(data$ego) - lastone #figure out how many are missing after that
    visitnw <- add_vertices(visitnw, numtoadd) #add these-- they'll be added as higher numbers
  }
  
  if(vcount(mealnw) < length(data$ego)){
    lastone <- max(V(mealnw)) #find the last vertex added
    numtoadd <- length(data$ego) - lastone #figure out how many are missing after that
    mealnw <- add_vertices(mealnw, numtoadd) #add these-- the'll be added as higher numbers
  }
  
  if(vcount(borrownw) < length(data$ego)){
    lastone <- max(V(borrownw)) #find the last vertex added
    numtoadd <- length(data$ego) - lastone #figure out how many are missing after that
    borrownw <- add_vertices(borrownw, numtoadd) #add these-- the'll be added as higher numbers
  }
  
  if(vcount(rumornw) < length(data$ego)){
    lastone <- max(V(rumornw)) #find the last vertex added
    numtoadd <- length(data$ego) - lastone #figure out how many are missing after that
    rumornw <- add_vertices(rumornw, numtoadd) #add these-- the'll be added as higher numbers
  }
  
  if(vcount(unionnw) < length(data$ego)){
    lastone <- max(V(unionnw)) #find the last vertex added
    numtoadd <- length(data$ego) - lastone #figure out how many are missing after that
    unionnw <- add_vertices(unionnw, numtoadd) #add these-- the'll be added as higher numbers
  }
  
  if(vcount(spokerefnw) < length(data$ego)){
    lastone <- max(V(spokerefnw)) #find the last vertex added
    numtoadd <- length(data$ego) - lastone #figure out how many are missing after that
    spokerefnw <- add_vertices(spokerefnw, numtoadd) #add these-- the'll be added as higher numbers
  }
  
  ## Add name and household number and household count to be node
  ## attributes in each network
  ## THIS IS THE PART THAT DIFFERS FROM OTHER CODE FILES
  ## We remove the hhnum attribute (it is sensitive to data order)
  ## We add an ego attribute
  ## And we assign the hhnum to the node's name
  

  V(visitnw)$hhcount <- hhcount
  V(visitnw)$ego <- data$ego
  V(visitnw)$name <- data$hhnum
  
  #V(mealnw)$hhnum <- 1:vcount(mealnw)
  #V(mealnw)$hhnum_check <- data$hhnum
  #V(mealnw)$name <- data$ego
  V(mealnw)$hhcount <- hhcount
  V(mealnw)$ego <- data$ego
  V(mealnw)$name <- data$hhnum
  
  V(borrownw)$hhcount <- hhcount
  V(borrownw)$ego <- data$ego
  V(borrownw)$name <- data$hhnum
  
  V(rumornw)$hhcount <- hhcount
  V(rumornw)$ego <- data$ego
  V(rumornw)$name <- data$hhnum
  
  V(unionnw)$hhcount <- hhcount
  V(unionnw)$ego <- data$ego
  V(unionnw)$name <- data$hhnum
  
  V(spokerefnw)$hhcount <- hhcount
  V(spokerefnw)$ego <- data$ego
  V(spokerefnw)$name <- data$hhnum
  
  ## Add other attributes in addatts function built above
  
  visitnw <- addatts(visitnw, data = data)
  mealnw <- addatts(mealnw, data = data)
  borrownw <- addatts(borrownw, data = data)
  rumornw <- addatts(rumornw, data = data)
  unionnw <- addatts(unionnw, data = data)
  spokerefnw <- addatts(spokerefnw, data = data)
  
  out <- list("visit" = visitnw, 
              "meal" = mealnw, 
              "borrow" = borrownw, 
              "rumor" = rumornw, 
              "union" = unionnw, 
              "spokeref" = spokerefnw) #collect all the networks for a village in a list
  
  return(out) #return the list of all of a village's networks
  
}


## Now use the function above to build all the networks for each
## of the 12 villages.

## Each village's networks will be stored in an object called v#n.  
## This is a list of village #'s seven networks (visit, meal, 
## borrow, rumor, union, willtell, and spokeref)

v1n <- makehhnws(data = v1, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v2n <- makehhnws(data = v2, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v3n <- makehhnws(data = v3, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v4n <- makehhnws(data = v4, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v5n <- makehhnws(data = v5, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v6n <- makehhnws(data = v6, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v7n <- makehhnws(data = v7, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v8n <- makehhnws(data = v8, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v9n <- makehhnws(data = v9, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v10n <- makehhnws(data = v10, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v11n <- makehhnws(data = v11, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)
v12n <- makehhnws(data = v12, directed = TRUE, rmduplinks = TRUE, rmloops = TRUE)

## We can grab any one of the networks with the $, for instance v5n$visit grabs
## the visit layer in village 5.  

## That will get cumbersome; it will be convenient to have the network we'll use
## frequently, the spoke refugees network, re-labeled something shorter
## than this.  So, pull out the spokeref networks and rename sr# for convenience:

sr1 <- v1n$spokeref
sr2 <- v2n$spokeref
sr3 <- v3n$spokeref
sr4 <- v4n$spokeref
sr5 <- v5n$spokeref
sr6 <- v6n$spokeref
sr7 <- v7n$spokeref
sr8 <- v8n$spokeref
sr9 <- v9n$spokeref
sr10 <- v10n$spokeref
sr11 <- v11n$spokeref
sr12 <- v12n$spokeref

## Let's do the same for the union networks.
## Pull out the union networks and rename g# for convenience:

g1 <- v1n$union
g2 <- v2n$union
g3 <- v3n$union
g4 <- v4n$union
g5 <- v5n$union
g6 <- v6n$union
g7 <- v7n$union
g8 <- v8n$union
g9 <- v9n$union
g10 <- v10n$union
g11 <- v11n$union
g12 <- v12n$union

## Now all the networks are created, and the spokeref and unions are renamed.
## The node names are the household numbers.

V(g3)[1] #household 102
V(g3)$ego[1] #100119 is the ego name of the 1st node in v3's network

which(D$village_number == 3 & D$ego == "100119") #the 1st row of the data
D$hhnum[1] #yep, household 102






###############
##
## Make Table 1: Aggregated social network by village

## We can streamline the calculation of all the network statistics we want
## for all 12 village's union networks with a function.

## So, let's build a function that calculates a bunch of network statistics and puts
## them in a matrix for any list of network objects that we feed it.

nwstats1 <- function(nwlist){ #take a list of networks and build a table of network stats
  outmat <- matrix(data = NA, nrow = length(nwlist), ncol = 11) #add a row for each network 
  
  for(i in 1:length(nwlist)){
    nw <- nwlist[[i]]
    outmat[i,1] <- vcount(nw) #number of nodes
    outmat[i,2] <- ecount(nw) #number of links
    outmat[i,3] <- mean(degree(nw)) #average total degree
    outmat[i,4] <- sort(degree(nw, mode = "out"), decreasing = TRUE)[1] #largest out-degree
    outmat[i,5] <- sort(degree(nw, mode = "in"), decreasing = TRUE)[1] #largest in-degree
    outmat[i,6] <- sum(degree(nw) == 0) #total number with degree 0
    outmat[i,7] <- sum(degree(nw, mode = "out") == 0) #total number with out-degree 0
    outmat[i,8] <- sum(degree(nw, mode = "in") == 0) #total number with in-degree 0
    outmat[i,9] <- transitivity(nw) #average transitivity
    outmat[i,10] <- components(nw)$no #number of components
    outmat[i,11] <- sort(components(nw)$csize, decreasing=TRUE)[1]/vcount(nw) #size of the largest component
  }
  colnames(outmat) <- c("Nodes", "Links", "Mean Deg", "Max OD", "Max ID", "Isol", 
                        "0 OD", "0 ID", "Trans", "#Comp", "LgComp")
  return(outmat)
}

## Now use the function to make the matrix of network stats for the union network 
## in all 12 villages and call it out1

out1 <- nwstats1(list(g1, g2, g3, g4, g5, g6, 
                      g7, g8, g9, g10, g11, g12))

## Now generate the latex code to build this table
xtable(out1, caption = "Collapsed Multi-Layer Network")


##################
##
## Make Table 2: Four layers separated by village

## Ok now build this same table, but this time separate out the four layers
## and average across 12 villages.

## A new function will make this easier.
## The function again takes a list of networks and returns a matrix of 
## network stats:


nwstats1b <- function(nwlist){
  outmat <- matrix(data = NA, nrow = length(nwlist), ncol = 8)
  
  for(i in 1:length(nwlist)){
    nw <- nwlist[[i]]
    outmat[i,1] <- vcount(nw)
    outmat[i,2] <- ecount(nw)
    outmat[i,3] <- mean(degree(nw))
    outmat[i,4] <- sort(degree(nw, mode = "in"), decreasing = TRUE)[1]
    outmat[i,5] <- sum(degree(nw, mode = "out") == 0)
    outmat[i,6] <- sum(degree(nw, mode = "in") == 0)
    outmat[i,7] <- transitivity(nw)
    outmat[i,8] <- sort(components(nw)$csize, decreasing = TRUE)[1]/vcount(nw)
  }
  colnames(outmat) <- c("Nodes", "Links", "Mean Deg", "Max ID", 
                        "0 OD", "0 ID", "Trans", "LgComp")
  return(outmat)
}

## Now calculate these stats for the layers separately:

outmeal <- nwstats1b(list(v1n$meal, v2n$meal, v3n$meal, v4n$meal, v5n$meal, v6n$meal,
                          v7n$meal, v8n$meal, v9n$meal, v10n$meal, v11n$meal, v12n$meal))

outvisit <- nwstats1b(list(v1n$visit, v2n$visit, v3n$visit, v4n$visit, v5n$visit, v6n$visit,
                           v7n$visit, v8n$visit, v9n$visit, v10n$visit, v11n$visit, v12n$visit))

outrumor <- nwstats1b(list(v1n$rumor, v2n$rumor, v3n$rumor, v4n$rumor, v5n$rumor, v6n$rumor,
                           v7n$rumor, v8n$rumor, v9n$rumor, v10n$rumor, v11n$rumor, v12n$rumor))

outborrow <- nwstats1b(list(v1n$borrow, v2n$borrow, v3n$borrow, v4n$borrow, v5n$borrow, v6n$borrow,
                            v7n$borrow, v8n$borrow, v9n$borrow, v10n$borrow, v11n$borrow, v12n$borrow))


##Average these across the 12 villages and collect the averages in a matrix called layers
layers <- rbind(apply(outmeal, 2, mean),
                apply(outvisit, 2, mean),
                apply(outrumor, 2, mean),
                apply(outborrow, 2, mean))

rownames(layers) <- c("Meal", "Visit", "Rumor", "Borrow") #label the rows

## We would like some of these to be rounded to zero decimal places, but others
## to have 2 decimals.  We can do that by building a matrix that will tell 
## x table how to round each entry in the matrix.  
## Call tht matrix diglayers and initialize it to have 0s everywhere.

diglayers <- matrix(data = 0, nrow = 4, ncol = 9) #set up the digits for the xtable

## Now replace the 0 with 2 for columns 4, 8, and 9, the ones we want to display
## with 2 instead of 0 digits

diglayers[,c(4,8,9)] <- 2

## Now generate the latex code for the table using the matrix of digits

xtable(layers, caption = "Characteristics of each of the four layers averaged over the 12 villages",
       digits = diglayers)

## Tidy our workspace
rm(diglayers, outmeal, outvisit, outrumor, outborrow)





#########
## --------------- OLS Models -------------------------------------------------
## Appendix Table A2, Treatment Effect on Network Link Usage

## First the simple ATE model, regressing the number of discussion partners
## on treatment status

model_basic_ate <- lm(D$numspokeref ~ D$trt, data = D)
summary(model_basic_ate)

## Now create outcome variables that count the number of discussion partners
## who also appear as neighbors in each network layer

D$numbervisits <- 0
D$numbermeals <- 0
D$numberborrows <- 0
D$numberrumors <- 0


# Loop through rows
for (i in 1:nrow(D)) {
  # Count occurrences for numbervisits
  D$numbervisits[i] <- sum(!is.na(D[i, paste0("visit", 1:5)]) & D[i, paste0("visit", 1:5)] %in% D[i, paste0("spokeref", 1:5)])
  
  # Count occurrences for numbermeals
  D$numbermeals[i] <- sum(!is.na(D[i, paste0("meal", 1:5)]) & D[i, paste0("meal", 1:5)] %in% D[i, paste0("spokeref", 1:5)])
  
  # Count occurrences for numberborrows
  D$numberborrows[i] <- sum(!is.na(D[i, paste0("borrow", 1:5)]) & D[i, paste0("borrow", 1:5)] %in% D[i, paste0("spokeref", 1:5)])
  
  # Count occurrences for numberrumors
  D$numberrumors[i] <- sum(!is.na(D[i, paste0("rumor", 1:5)]) & D[i, paste0("rumor", 1:5)] %in% D[i, paste0("spokeref", 1:5)])
}

## Regress each of these DVs on treatment status

numbervisitsols <- lm(D$numbervisits ~ D$trt, data = D)
summary(numbervisitsols) ## Treated more likely to discuss refugees with more visit links. 

numbermealsols <- lm(D$numbermeals ~ D$trt, data = D)
summary(numbermealsols) ## Treated more likely to discuss refugees with more meal links. 

numberborrowsols <- lm(D$numberborrows ~ D$trt, data = D)
summary(numberborrowsols) ## Treated more likely to discuss refugees with more borrow links. 

numberrumorsols <- lm(D$numberrumors ~ D$trt, data = D)
summary(numberrumorsols) ## Null for rumors, I think that was expected based on previous conversations? 

## And build Appendix Table A2

model_list <- list(model_basic_ate, numbermealsols, numbervisitsols, 
                   numberrumorsols, numberborrowsols)
stargazer(model_list)


#####
## Now repeat, including village fixed effects here 
## (since we've shown there's lots of variation at the village level)


## Adding village fixed effects to all of the above:
model_basic_ate_vfe <- lm(numspokeref ~ D$trt + as.character(D$village_number), data = D)
summary(model_basic_ate_vfe)


numbervisitsols_vfe <- lm(D$numbervisits ~ D$trt + as.character(D$village_number), data = D)

numbermealsols_vfe <- lm(D$numbermeals ~ D$trt + as.character(D$village_number), data = D)

numberborrowsols_vfe <- lm(D$numberborrows ~ D$trt + as.character(D$village_number), data = D)


numberrumorsols_vfe <- lm(D$numberrumors ~ D$trt + as.character(D$village_number), data = D)


model_list_vfe <- list(model_basic_ate_vfe, numbermealsols_vfe, numbervisitsols_vfe, 
                       numberrumorsols_vfe, numberborrowsols_vfe)


# And build Paper Table 4
stargazer(model_list_vfe,
          keep.stat = c("n", "adj.rsq"),
          omit = c("village_number"),
          omit.labels = c("Village FE"))





###################
## Next, let's implement a permutation test to verify that the difference in means
## is robust to network dependencies

## Run the QAP, which here is equivalent to randomly reassigning T and C
## within each village (hold the matrix of sender T v. C constant while shuffling
## the rows and columns in the who listed whom as conversation partners adjecency
## matrix, which is then used to compute out-degree)

## Let's permute the treatment assignment 1000 times within each village and
## record the difference in means between T and C with respect to each of the 
## five outcomes: numspokeref, numbermeals, numbervisits, numberrumors, numberborrows


nsims <- 1000
set.seed(20240415)

numvlgs <- length(unique(D$village_number))

sim_reg_coefs <- matrix(data = NA, nrow = nsims, ncol = 10)
colnames(sim_reg_coefs) <- c("all", "visit", "meal", "borrow", "rumor",
                             "all_fe", "visit_fe", "meal_fe", "borrow_fe", "rumor_fe")

D_sm<- D %>% select(village_number, trt, numspokeref, numbervisits, numbermeals,
                    numberborrows, numberrumors)


for (i in 1:nsims){
  
  ## Make an empty vector the length of our dataset
  simtrt <- rep(NA, nrow(D_sm))
  
  ## For each village, randomly reassign the treatment vector
  
  for(j in 1:numvlgs){
    
    ## Identify the rows that are part of village j
    thevlg <- which(D_sm$village_number == j)
    
    ## Grab the treatment values for this village
    thetrt <- D_sm$trt[thevlg]
    
    ## Randomly shuffle these values
    newtrt <- sample(thetrt, length(thetrt), replace = FALSE)
    
    ## Assign the shuffled values to the right spots in the simtrt vector
    simtrt[thevlg] <- newtrt
  }
  
  ## Now simtrt is a simulated treatment vector the length of the D_sm dataset 
  ## Add it to our dataset
  
  D_sm$simtrt <- simtrt
  
  ## Regress our DVs on this one and store the coefficient on simtrt (in spot 2, intercept is 1)
  sim_reg_coefs[i,1] <- coef(lm(D_sm$numspokeref ~ D_sm$simtrt, data = D_sm))[2]
  sim_reg_coefs[i,2] <- coef(lm(D_sm$numbervisits ~ D_sm$simtrt, data = D_sm))[2]
  sim_reg_coefs[i,3] <- coef(lm(D_sm$numbermeals ~ D_sm$simtrt, data = D_sm))[2]
  sim_reg_coefs[i,4] <- coef(lm(D_sm$numberborrows ~ D_sm$simtrt, data = D_sm))[2]
  sim_reg_coefs[i,5] <- coef(lm(D_sm$numberrumors ~ D_sm$simtrt, data = D_sm))[2]
  
  ## And repeat, adding village fixed effects:
  sim_reg_coefs[i,6] <- coef(lm(D_sm$numspokeref ~ D_sm$simtrt + as.character(D$village_number), data = D_sm))[2]
  sim_reg_coefs[i,7] <- coef(lm(D_sm$numbervisits ~ D_sm$simtrt + as.character(D$village_number), data = D_sm))[2]
  sim_reg_coefs[i,8] <- coef(lm(D_sm$numbermeals ~ D_sm$simtrt + as.character(D$village_number), data = D_sm))[2]
  sim_reg_coefs[i,9] <- coef(lm(D_sm$numberborrows ~ D_sm$simtrt + as.character(D$village_number), data = D_sm))[2]
  sim_reg_coefs[i,10] <- coef(lm(D_sm$numberrumors ~ D_sm$simtrt + as.character(D$village_number), data = D_sm))[2]
  
  if(i %% 100 == 0) print(i)
    
}
  
## And now we compare our output to the observed regression coefficients

coef_num <- coef(model_basic_ate)[2]
coef_visits <- coef(numbervisitsols)[2]
coef_meals <- coef(numbermealsols)[2]
coef_borrows <- coef(numberborrowsols)[2]
coef_rumors <- coef(numberrumorsols)[2]

coef_num_fe <- coef(model_basic_ate_vfe)[2]
coef_visits_fe <- coef(numbervisitsols_vfe)[2]
coef_meals_fe <- coef(numbermealsols_vfe)[2]
coef_borrows_fe <- coef(numberborrowsols_vfe)[2]
coef_rumors_fe <- coef(numberrumorsols_vfe)[2]


## And now assess significance
sum(sim_reg_coefs[,1] >= coef_num) #0
sum(sim_reg_coefs[,2] >= coef_visits) #7
sum(sim_reg_coefs[,3] >= coef_meals) #1
sum(sim_reg_coefs[,4] >= coef_borrows) #1
sum(sim_reg_coefs[,5] >= coef_rumors) #547

sum(sim_reg_coefs[,6] >= coef_num_fe) #0
sum(sim_reg_coefs[,7] >= coef_visits_fe) #7
sum(sim_reg_coefs[,8] >= coef_meals_fe) #1
sum(sim_reg_coefs[,9] >= coef_borrows_fe) #1
sum(sim_reg_coefs[,10] >= coef_rumors_fe) #545
## FE models very similar, which we would expect since we perform the 
## randomization within villages; they should be independent across


## And calculate our two-sided QAP scores 
sum(abs(sim_reg_coefs[,1]) >= abs(coef_num))/1000 #.000
sum(abs(sim_reg_coefs[,2]) >= abs(coef_visits))/1000 #.017
sum(abs(sim_reg_coefs[,3]) >= abs(coef_meals))/1000 #.001
sum(abs(sim_reg_coefs[,4]) >= abs(coef_borrows))/1000 #.004
sum(abs(sim_reg_coefs[,5]) >= abs(coef_rumors))/1000 #.833

sum(abs(sim_reg_coefs[,6]) >= abs(coef_num_fe))/1000 #.000
sum(abs(sim_reg_coefs[,7]) >= abs(coef_visits_fe))/1000 #.012
sum(abs(sim_reg_coefs[,8]) >= abs(coef_meals_fe))/1000 #.001
sum(abs(sim_reg_coefs[,9]) >= abs(coef_borrows_fe))/1000 #.003
sum(abs(sim_reg_coefs[,10]) >= abs(coef_rumors_fe))/1000 #.879



sim_reg_coefs <- as.data.frame(sim_reg_coefs)


## Let's plot the distributions one at a time

sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(all), fill = "gold1", alpha = .8)+
  geom_vline(xintercept = coef_num, linetype = "dashed")+
  labs(title = "QAP, All Conversation Partners", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_all.pdf, 6x6

sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(all_fe), fill = "gold1", alpha = .8)+
  geom_vline(xintercept = coef_num_fe, linetype = "dashed")+
  labs(title = "QAP, All Conversation Partners, FE", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_all_fe.pdf, 6x6

sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(visit), fill = "lightskyblue", alpha = .8)+
  geom_vline(xintercept = coef_visits, linetype = "dashed")+
  labs(title = "QAP, Visit Layer", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_visit.pdf, 6x6

sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(visit_fe), fill = "lightskyblue", alpha = .8)+
  geom_vline(xintercept = coef_visits_fe, linetype = "dashed")+
  labs(title = "QAP, Visit Layer, FE", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_visit_fe.pdf

sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(meal), fill = "lightgreen", alpha = .8)+
  geom_vline(xintercept = coef_meals, linetype = "dashed")+
  labs(title = "QAP, Meal Layer", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_meal.pdf, 6x6

sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(meal_fe), fill = "lightgreen", alpha = .8)+
  geom_vline(xintercept = coef_meals_fe, linetype = "dashed")+
  labs(title = "QAP, Meal Layer, FE", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_meal_fe.pdf

sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(borrow), fill = "thistle3", alpha = .8)+
  geom_vline(xintercept = coef_borrows, linetype = "dashed")+
  labs(title = "QAP, Borrow Layer", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_borrow.pdf, 6x6

sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(borrow_fe), fill = "thistle3", alpha = .8)+
  geom_vline(xintercept = coef_borrows_fe, linetype = "dashed")+
  labs(title = "QAP, Borrow Layer, FE", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_borrow_fe.pdf


sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(rumor), fill = "lightsalmon", alpha = .8)+
  geom_vline(xintercept = coef_rumors, linetype = "dashed")+
  labs(title = "QAP, Rumor Layer", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_rumor.pdf, 6x6

sim_reg_coefs %>%
  ggplot()+
  geom_density(aes(rumor_fe), fill = "lightsalmon", alpha = .8)+
  geom_vline(xintercept = coef_rumors_fe, linetype = "dashed")+
  labs(title = "QAP, Rumor Layer, FE", 
       x = "Regression Coefficients from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_rumor_fe.pdf





#####
## Table 5
## Comparing number spoke to, t v. c by village:


numspoketvc <- matrix(data = NA, nrow = 12, ncol = 2)
colnames(numspoketvc) <- c("T", "C")

for(i in 1:length(unique(D$village_number))){ #villages range from 1 to 12
  thisround <- which(D$village_number == i)
  vl <- D[thisround,]
  
  thet <- which(vl$trt == 1)
  thec <- which(vl$trt == 0)
  
  numspoketvc[i,1] <- mean(vl$numspokeref[thet])
  numspoketvc[i,2] <- mean(vl$numspokeref[thec])
  
}


xtable(numspoketvc, digits = 1, caption = "Number of Discussion Partners by Treatment Condition")

## And the pooled total:

mean(D$numspokeref[D$trt == 1]) #2.1
mean(D$numspokeref[D$trt == 0]) #1.6







##########################
## Table 6

## Explore the overlap in networks
## We want to know how many names listed as refugee discussion partners 
## (links in the spokeref network) are also listed in each of the 
## layers of the social network measured earlier as well as in the union
## of the layers.

## To do this, first, collect the (deidentified) names relevant to each of the layers.
## We'll first make a vector of the right variable names in the dataset

mlnames <- c("meal1", "meal2", "meal3", "meal4", "meal5")
vsnames <- c("visit1", "visit2", "visit3", "visit4", "visit5")
rmnames <- c("rumor1", "rumor2", "rumor3", "rumor4", "rumor5")
bwnames <- c("borrow1", "borrow2", "borrow3", "borrow4", "borrow5")

dtnames <- c("spokeref1", "spokeref2", "spokeref3", "spokeref4", "spokeref5")

## And then we'll make dataframes of the names listed for each layer by referencing
## these variable names.  

ml <- D[,names(D) %in% mlnames]
vs <- D[,names(D) %in% vsnames]
rm <- D[,names(D) %in% rmnames]
bw <- D[,names(D) %in% bwnames]
dt <- D[,names(D) %in% dtnames]

## We can combine the four layers' dataframes to make the dataframe of names
## in the union network.

union <- cbind(ml, vs, rm, bw)

## Now let's write a function that takes two dataframes of names
## and returns the count and the proportion of names appearing in one
## row of dataframe a that area also in that row of dataframe b.
## This will return the overlap of two networks.
## (in otherwords, how much a overlaps with b)

ainb <- function(a,b){ #a and b are dfs of name variables,
  outmat <- matrix(data = NA, nrow = nrow(a), ncol = 2)
  colnames(outmat) <- c("numin", "propin")
  
  for(i in 1:nrow(a)){
    anoblanks <- a[i,][which(a[i,] != "")]
    bnoblanks <- b[i,][which(b[i,] != "")]
    outmat[i,1] <- sum(anoblanks %in% bnoblanks)
    outmat[i,2] <- sum(anoblanks %in% bnoblanks)/length(anoblanks)
    
  }
  return(outmat)
  
}

## To make sure it's working, let's look at the names in the spokeref network
## that are in the meal layer.

test <- ainb(dt, ml)
dim(test) #1 row for each of the 1604 observations, 2 columns.
head(test) #reports the number and proportion of names in dt that are in ml.

## Ok, now for each network, let's store the extent of overlap of each layer with didtell
## as a variable in the network.
## This is the number of names in the did tell (spokeref) network that also appear
## in that person's meal, visit, rumor, borrow, and union networks.

D$dtinmeal <- ainb(dt, ml)[,1] #1st column because we're saving numin and discarding propin
D$dtinvisit <- ainb(dt, vs)[,1]
D$dtinrumor <- ainb(dt, rm)[,1]
D$dtinborrow <- ainb(dt, bw)[,1]
D$dtinunion <- ainb(dt, union)[,1]


## So now, let's make a table that shows how the names people offered as the people
## with whom they spoke about refugees (spokeref) appear in the union network 
## and each of the layers.  This will be conditional on having told anyone 
## (the number of names offered for spokeref, numspokeref, is positive).
## This will make Paper Table 6.

overlapmat4 <- matrix(data = NA, nrow = 12, ncol = 5)
colnames(overlapmat4) <- c("Tot in NA", "#inMeal", "#inVisit", "#inRumor", "#inBorrow")

for(i in 1:length(unique(D$village_number))){ #villages range from 1 to 12
  thisround <- which(D$village_number == i)
  vl <- D[thisround,]
  didtells <- which(vl$numspokeref > 0)
  
  overlapmat4[i,1] <- mean(vl$dtinunion[didtells])
  overlapmat4[i,2] <- mean(vl$dtinmeal[didtells])
  overlapmat4[i,3] <- mean(vl$dtinvisit[didtells])
  overlapmat4[i,4] <- mean(vl$dtinrumor[didtells])
  overlapmat4[i,5] <- mean(vl$dtinborrow[didtells])
  
}


xtable(overlapmat4, digits = 2, caption = "Breakdown of discussion partners by layers of network")

## And the total averaged over the villages:
round(apply(overlapmat4, 2, mean),2)




###############
##
## Now we can build Appendix Table A1.
## In each village, breakdown of spokeref ties 
## across four villages.  How many people said they talked to anyone,
## what proportion of the village this is, and the proportion of the people
## listed as the people with whom they spoke about refugees also appear
## in that ego's union network 

spokerefmat1 <- matrix(data = NA, nrow = 12, ncol = 3)
colnames(spokerefmat1) <- c("#DT>0", "PropDT>0", "AnyInNW")

for(i in 1:length(unique(D$village_number))){ #villages range from 1 to 12
  thisround <- which(D$village_number == i)
  vl <- D[thisround,]
  didtells <- which(vl$numspokeref > 0)
  spokerefmat1[i,1] <- length(didtells)
  spokerefmat1[i,2] <- length(didtells)/nrow(vl)
  spokerefmat1[i,3] <- sum(vl$dtinunion[didtells] > 0)/length(didtells) #proportion of people who told at least 1 for whom at least 1 was in union
  
}


digspokerefmat1 <- matrix(data = 0, nrow = 12, ncol = 4) #set the digits to be 0 for all table entries
digspokerefmat1[,c(3,4)] <- 2 #make columns 3 and for be 2 digits


## Appendix Table A1:
xtable(spokerefmat1, caption = "Who Told Anyone, Were They Network Neighbors?", digits = digspokerefmat1)







######################################################
######################################################
## Now build a dataset organized by links



################################
## Link dataset

## The full dataset D and each of the 12 village ones v1...v12 are all
## organized by individual.  These are our egos in the network, and the
## alter information is listed as additional columns.  We used these to build
## edgelists, make networks, convert them to household ones, etc.

## Now we want to reorganize the data to have one row for every link in
## one of the networks (union/layers, spokeref).

## To do this, we will first back out our edgelists from the networks in 
## terms of household numbers.  We will add a column indicating the village number
## and another of 1s that will become the indicator for appearing in that network.
## We will then merge all of these edgelsits in terms of the household
## number of the ego, alter, and village number, keeping all rows.
## Then we will add ego attributes and then alter attributes from our individual
## level dataframe.

## Step 1: create datagrames from edgelists for each of meals, visits, borrow, rumors,
## union, and spokeref in all 12 villages.

## We can write a function that starts with a network, extracts the edgelist
## in terms of the household number (rather than the individual's name),
## adds a column indicating village number and adds a column of 1s
## with a variable name we give it.


## vilnum is a number (the village number corresponding to the network) 
## indicname is a string to name the indicator variable in terms of the netowrk supplied
makehheldf <- function(nw, vilnum, indicname){ #grabs el plus adds specific variables from larger dataset with attributes in it
  
  #grab edgelist in terms of names, which in the network creation above are
  #household numbers in the data unique to each village
  getel <- get.edgelist(nw, names = TRUE) 
  
  #add a colummn of 1s
  elplusindic <- cbind(getel, 1) 
  
  #add names for these columns using the name provided in indicname for the 3rd
  colnames(elplusindic) <- c("ego_hh", "alter_hh", indicname) 
  
  #turn this matrix into a data frame
  newdata <- as.data.frame(elplusindic) 
  
  #add a villag_number variable that enters the number provided in vilnum
  newdata$village_number <- vilnum 
  
  return(newdata)
}


## Let's make sure it's working
test <- makehheldf(g1, vilnum = 1, indicname = "in_union")
dim(test) #799x4
head(test) #good, all 4 expected columns.  And is 799 right?
ecount(g1) #yep, there are 799 links in village 1's union network g1

## One more test:
test2 <- makehheldf(v2n$meal, vilnum = 2, indicname = "in_meal")
dim(test2) #241x4
head(test2) #good, village number is showing up 2 as expected
ecount(v2n$meal) #good, there are 241 links in village 2's meal network

rm(test, test2)


## Now we need to use this function to make these dataframes for
## all  6 networks (meal, visit, rumor, borrow, union, spokeref) for 
## all 12 villages.

## Let's make all of meal for all 12 villages:

numvlgs <- 12

mealnwlist <- list(v1n$meal, v2n$meal, v3n$meal, v4n$meal,
                   v5n$meal, v6n$meal, v7n$meal, v8n$meal,
                   v9n$meal, v10n$meal, v11n$meal, v12n$meal)

## Note, this loop assumes the villages are in numeric order in the network list.  
## When i is 3, it assumes it's working on Village 3.

## Get started by making the dataframe for the first village
meallinks <- makehheldf(mealnwlist[[1]], vilnum = 1, indicname = "in_meal") #

for(i in 2:numvlgs){ #for each of villages 2- 12,
  
  thenw <- mealnwlist[[i]] #take that village's network
  newdata <- makehheldf(thenw, vilnum = i, indicname = "in_meal") #make the village dataframe
  
  ## Now add the rows to the bottom of the current version of meallinks dataframe
  meallinks <- rbind(meallinks, newdata)
  
}

## Now meallinks is a dataframe with all 12 villages' meallinks in it

table(meallinks$village_number)

## And we can make sure these counts are right by quickly using lapply
## to apply the function that counts edges in a network, e count,
## to every item in our meal network list.  Unlist will report the results 
## as a vector (instead of a default list).  

unlist(lapply(mealnwlist, ecount))

## Good, all meal links in each village are accounted for.

rm(mealnwlist)

## Let's do the same with the visit layer:

visitnwlist <- list(v1n$visit, v2n$visit, v3n$visit, v4n$visit,
                    v5n$visit, v6n$visit, v7n$visit, v8n$visit,
                    v9n$visit, v10n$visit, v11n$visit, v12n$visit)

## Start with the visit network for village 1
visitlinks <- makehheldf(visitnwlist[[1]], vilnum = 1, indicname = "in_visit") #note we changed the label to in_visit

for(i in 2:numvlgs){ # repeat for each of villages 2- 12,
  
  thenw <- visitnwlist[[i]] #take that village's network
  newdata <- makehheldf(thenw, vilnum = i, indicname = "in_visit") #make the village dataframe
  
  ## Now add the rows to the bottom of the current version of dataframe
  visitlinks <- rbind(visitlinks, newdata)
  
}

## Spot check the counts:
table(visitlinks$village_number)
unlist(lapply(visitnwlist, ecount)) #good

rm(visitnwlist)

## And repeating with the rumor layer


rumornwlist <- list(v1n$rumor, v2n$rumor, v3n$rumor, v4n$rumor,
                    v5n$rumor, v6n$rumor, v7n$rumor, v8n$rumor,
                    v9n$rumor, v10n$rumor, v11n$rumor, v12n$rumor)

## Start with the rumor network for village 1
rumorlinks <- makehheldf(rumornwlist[[1]], vilnum = 1, indicname = "in_rumor") 

for(i in 2:numvlgs){ # repeat for each of villages 2- 12,
  
  thenw <- rumornwlist[[i]] #take that village's network
  newdata <- makehheldf(thenw, vilnum = i, indicname = "in_rumor") #make the village dataframe
  
  ## Now add the rows to the bottom of the current version of  dataframe
  rumorlinks <- rbind(rumorlinks, newdata)
  
}

table(rumorlinks$village_number)
unlist(lapply(rumornwlist, ecount)) #good

rm(rumornwlist)

## And the borrow layer:

borrownwlist <- list(v1n$borrow, v2n$borrow, v3n$borrow, v4n$borrow,
                     v5n$borrow, v6n$borrow, v7n$borrow, v8n$borrow,
                     v9n$borrow, v10n$borrow, v11n$borrow, v12n$borrow)

## Start with the rumor network for village 1
borrowlinks <- makehheldf(borrownwlist[[1]], vilnum = 1, indicname = "in_borrow") 

for(i in 2:numvlgs){ # repeat for each of villages 2- 12,
  
  thenw <- borrownwlist[[i]] #take that village's network
  newdata <- makehheldf(thenw, vilnum = i, indicname = "in_borrow") #make the village dataframe
  
  ## Now add the rows to the bottom of the current version of dataframe
  borrowlinks <- rbind(borrowlinks, newdata)
  
}

table(borrowlinks$village_number)
unlist(lapply(borrownwlist, ecount)) #good

rm(borrownwlist)

## And the union of all the layers (although this will pick up redundant info,
## merge will take care of it, and this is the easiest way to have an in_union
## indicator in the ultimate link dataset).


unionnwlist <- list(g1, g2, g3, g4,
                    g5, g6, g7, g8,
                    g9, g10, g11, g12)

## Start with the union network for village 1
unionlinks <- makehheldf(unionnwlist[[1]], vilnum = 1, indicname = "in_union") 

for(i in 2:numvlgs){ # repeat for each of villages 2- 12,
  
  thenw <- unionnwlist[[i]] #take that village's network
  newdata <- makehheldf(thenw, vilnum = i, indicname = "in_union") #make the village dataframe
  
  ## Now add the rows to the bottom of the current version of dataframe
  unionlinks <- rbind(unionlinks, newdata)
  
}

table(unionlinks$village_number)
unlist(lapply(unionnwlist, ecount)) #good


rm(unionnwlist)

## And finally, the spokeref networks:

spokerefnwlist <- list(sr1, sr2, sr3, sr4,
                       sr5, sr6, sr7, sr8,
                       sr9, sr10, sr11, sr12)

## Start with the spokeref network for village 1
spokereflinks <- makehheldf(spokerefnwlist[[1]], vilnum = 1, indicname = "in_spokeref") 

for(i in 2:numvlgs){ # repeat for each of villages 2- 12,
  
  thenw <- spokerefnwlist[[i]] #take that village's network
  newdata <- makehheldf(thenw, vilnum = i, indicname = "in_spokeref") #make the village dataframe
  
  ## Now add the rows to the bottom of the current version of dataframe
  spokereflinks <- rbind(spokereflinks, newdata)
  
}

table(spokereflinks$village_number)
unlist(lapply(spokerefnwlist, ecount)) #good

rm(spokerefnwlist)

## Ok, now we have 6 dataframes of links.   We can merge these by ego hh, alter hh, and village_number.

## Do this pair-wise.  First merge meal and visit link dataframes.
Links <- merge(meallinks, visitlinks, by = c("ego_hh", "alter_hh", "village_number"),
               all.x = TRUE, all.y=TRUE)

## Next merge this with borrow dataframe.
Links <- merge(Links, borrowlinks, by = c("ego_hh", "alter_hh", "village_number"),
               all.x = TRUE, all.y=TRUE)

## Next merge this with rumor dataframe.
Links <- merge(Links, rumorlinks, by = c("ego_hh", "alter_hh", "village_number"),
               all.x = TRUE, all.y=TRUE)

dim(Links) #7805 x 7 so far. Now merge with union; this shouldn't add any unique rows
## because all links in the union necessarily are in at least one layer.

Links <- merge(Links, unionlinks, by = c("ego_hh", "alter_hh", "village_number"),
               all.x = TRUE, all.y=TRUE)

dim(Links) #7805 x 8.  Good.  Added a column but no new rows.

## And finally, merge with spokeref.

Links <- merge(Links, spokereflinks, by = c("ego_hh", "alter_hh", "village_number"),
               all.x = TRUE, all.y=TRUE)

dim(Links) #8746 x 9.  
head(Links) #Excellent.  Now we need to replace the NAs in the in_XX variables
## with 0s, because these mean the link is not in that layer.

Links$in_meal[which(is.na(Links$in_meal))] <- 0
Links$in_visit[which(is.na(Links$in_visit))] <- 0
Links$in_borrow[which(is.na(Links$in_borrow))] <- 0
Links$in_rumor[which(is.na(Links$in_rumor))] <- 0
Links$in_union[which(is.na(Links$in_union))] <- 0
Links$in_spokeref[which(is.na(Links$in_spokeref))] <- 0

## And tidy the workspace

rm(meallinks, visitlinks, borrowlinks, rumorlinks, unionlinks, spokereflinks)



## Next we need to add some ego attributes, for example:
atts <- c("been_refugee",
          "refugee_meet",
          "threat1", 
          "trt",
          "farmer",
          "resp_language",
          "religion",
          "important_surr_refug")


## Now we need to create a dataset that includes the household number,
## village number, and these attributes.  



## First grab these attributes in their own dataset
egoatts <- D[,which(names(D) %in% atts)]

## Now add "_ego" to the names of these variables
names(egoatts) <- paste(names(egoatts), "ego", sep = "_")

## And add the household number which we'll call ego_hh
egoatts$ego_hh <- D$hhnum
## And the village number
egoatts$village_number <- D$village_number

## Now merge into Links dataset by ego hh and village number
Links <- merge(Links, egoatts, by = c("ego_hh", "village_number"), all.x = TRUE, all.y = FALSE)

rm(egoatts)

## Repeat this same set of steps, this time calling the household number alter_hh
## adding "_alter" to the variable names, and merging by alter_hh and village number

alteratts <- D[,which(names(D) %in% atts)]

## Now add "_ego" to the names of these variables
names(alteratts) <- paste(names(alteratts), "alter", sep = "_")

## And add the household number which we'll call alter_hh
alteratts$alter_hh <- D$hhnum
## And the village number
alteratts$village_number <- D$village_number

## Now merge into Links dataset by ego hh and village number
Links <- merge(Links, alteratts, by = c("alter_hh", "village_number"), all.x = TRUE, all.y = FALSE)

dim(Links) #8746 x 25
names(Links)

rm(alteratts)


#######
## Now Links is the dataframe where each row is a link


#####
## Now we can make Paper Table 7

## Let's compare links in the union social network that were and were
## not used to discuss refugees (that is, that were and were not in spokeref)


## Make a matrix to hold the comparisons
srmat <- matrix(data = NA, nrow = 3, ncol = 11)

## Make an index to identify the rows in Links that are in the spokeref network
## and in the union social network
srs <- which(Links$in_spokeref == 1 & Links$in_union == 1)

## Make an index to identify the rows in Links that are NOT in the spokeref network
## but ARE in the union social network
notsrs <- which(Links$in_spokeref == 0 & Links$in_union == 1)

## We can now fill out each column of our 3 by 11 matrix with the comparisons 
## want:

srmat[1,1] <- length(srs) #number of links in union nw used to discuss refs
srmat[2,1] <- length(notsrs) #number of links in unio nw not used to discuss refs

## Are they seeking out someone with relevant experience?
#### Has been refugee?
srmat[1,2] <- mean(Links$been_refugee_alter[srs], na.rm=TRUE) 
srmat[2,2] <- mean(Links$been_refugee_alter[notsrs], na.rm=TRUE)
srmat[3,2] <- t.test(Links$been_refugee_alter[srs], Links$been_refugee_alter[notsrs])$p.value


#### Knows refugee?
srmat[1,3] <- mean(Links$refugee_meet_alter[srs], na.rm=T)
srmat[2,3] <- mean(Links$refugee_meet_alter[notsrs], na.rm=T)
srmat[3,3] <- t.test(Links$refugee_meet_alter[srs], Links$refugee_meet_alter[notsrs])$p.value


#### Farmer?
srmat[1,4] <- mean(Links$farmer_alter[srs], na.rm=T)
srmat[2,4] <- mean(Links$farmer_alter[notsrs], na.rm=T)
srmat[3,4] <- t.test(Links$farmer_alter[srs], Links$farmer_alter[notsrs])$p.value

#### Feels a certain way about the topic?
srmat[1,5] <- mean(Links$threat1_alter[srs], na.rm=T)
srmat[2,5] <- mean(Links$threat1_alter[notsrs], na.rm=T)
srmat[3,5] <- t.test(Links$threat1_alter[srs], Links$threat1_alter[notsrs])$p.value #not sig

#### Finds the topic to be important?

srmat[1,6] <- mean(Links$important_surr_refug_alter[srs], na.rm=T)
srmat[2,6] <- mean(Links$important_surr_refug_alter[notsrs], na.rm=T)
srmat[3,6] <- t.test(Links$important_surr_refug_alter[srs], Links$important_surr_refug_alter[notsrs])$p.value #ok, importance of alter matters

## Is culturally similar to the respondent in terms of religion?

srmat[1,7] <- mean(Links$religion_ego[srs] == Links$religion_alter[srs], na.rm=T)
srmat[2,7] <- mean(Links$religion_ego[notsrs] == Links$religion_alter[notsrs], na.rm=T)
srmat[3,7] <- t.test(Links$religion_ego[srs] == Links$religion_alter[srs],
                     Links$religion_ego[notsrs] == Links$religion_alter[notsrs])$p.value #not sig, close-ish

#### In terms of language?

srmat[1,8] <- mean(Links$resp_language_ego[srs] == Links$resp_language_alter[srs], na.rm=T)
srmat[2,8] <- mean(Links$resp_language_ego[notsrs] == Links$resp_language_alter[notsrs], na.rm=T)
srmat[3,8] <- t.test(Links$resp_language_ego[srs] == Links$resp_language_alter[srs],
                     Links$resp_language_ego[notsrs] == Links$resp_language_alter[notsrs])$p.value

#### In terms of shared refugee experience?

srmat[1,9] <- mean(Links$been_refugee_ego[srs] == Links$been_refugee_alter[srs], na.rm=T)
srmat[2,9] <- mean(Links$been_refugee_ego[notsrs] == Links$been_refugee_alter[notsrs], na.rm=T)
srmat[3,9] <- t.test(Links$been_refugee_ego[srs] == Links$been_refugee_alter[srs],
                     Links$been_refugee_ego[notsrs] == Links$been_refugee_alter[notsrs])$p.value #not homoph w.r.t. ref status

#### In terms of shared views on the topic?
srmat[1,10] <- mean(Links$threat1_ego[srs] == Links$threat1_alter[srs], na.rm=T)
srmat[2,10] <- mean(Links$threat1_ego[notsrs] == Links$threat1_alter[notsrs], na.rm=T)
srmat[3,10] <- t.test(Links$threat1_ego[srs] == Links$threat1_alter[srs], Links$threat1_ego[notsrs] == Links$threat1_alter[notsrs])$p.value #sig at .1, seek out similarit

#### In terms of shared perception of importance?
srmat[1,11] <- mean(Links$important_surr_refug_ego[srs] == Links$important_surr_refug_alter[srs], na.rm=T)
srmat[2,11] <- mean(Links$important_surr_refug_ego[notsrs] == Links$important_surr_refug_alter[notsrs], na.rm=T)
srmat[3,11] <- t.test(Links$important_surr_refug_ego[srs] == Links$important_surr_refug_alter[srs],
                      Links$important_surr_refug_ego[notsrs] == Links$important_surr_refug_alter[notsrs]  )$p.value #sig at .5


colnames(srmat) <- c("Links", "Alter was ref", "Alter knows ref", "Alter farmer", "Alter views", "Alter interest", "Relig homoph", "Language homoph", 
                     "ref homoph", "views homoph", "interest homoph")

rownames(srmat) <- c("Link used to discuss", "Link not used to discuss", "p-value")


## And build Paper Table 6, first 10 rows

xtable(t(srmat), digits = 2)


###############
## Next, remake the standard errors using QAP for each of these comparisons

## We can simplify these by calculating the difference in means as simple ols regression
## coefficients. We'll need a dataframe that includes only links that are in the social 
## network, and a DV that takes 1 if that link was also listed as a spoke about refugees link

L <- Links %>% 
  filter(in_union == 1)

## Now L is the dataframe of Links that only includes the ones in at least one of the four layers
## and in_spokeref is the DV

table(L$in_spokeref)

## Now our p-values listed in Table 6 are equivalent to the ones from these regressions:
summary(lm(in_spokeref ~ been_refugee_alter, data = L))
summary(lm(in_spokeref ~ refugee_meet_alter, data = L))
summary(lm(in_spokeref ~ been_refugee_alter, data = L))
summary(lm(in_spokeref ~ threat1_alter, data = L))
summary(lm(in_spokeref ~ important_surr_refug_alter, data = L))

## Good, these replicate.  And now for the homophily ones, we need indicators of
## homophily

L$relig_homoph <- ifelse(L$religion_ego == L$religion_alter, 1, 0)
L$lang_homoph <- ifelse(L$resp_language_ego == L$resp_language_alter, 1, 0)
L$ref_homoph <- ifelse(L$been_refugee_ego == L$been_refugee_alter, 1, 0)
L$views_homoph <- ifelse(L$threat1_ego == L$threat1_alter, 1, 0)
L$interest_homoph <- ifelse(L$important_surr_refug_ego == L$important_surr_refug_alter, 1, 0)

summary(lm(in_spokeref ~ relig_homoph, data = L))
summary(lm(in_spokeref ~ lang_homoph, data = L))
summary(lm(in_spokeref ~ ref_homoph, data = L))
summary(lm(in_spokeref ~ views_homoph, data = L))
summary(lm(in_spokeref ~ interest_homoph, data = L))


## Next, let's add a more nuanced view of interest

## Let's make factor variables for homophily on different values

L$both_5 <- ifelse(L$views_homoph == 1 & L$threat1_ego == 5, 1, 0)
L$both_4 <- ifelse(L$views_homoph == 1 & L$threat1_ego == 4, 1, 0)
L$both_3 <- ifelse(L$views_homoph == 1 & L$threat1_ego == 3, 1, 0)
L$both_2 <- ifelse(L$views_homoph == 1 & L$threat1_ego == 2, 1, 0)
L$both_1 <- ifelse(L$views_homoph == 1 & L$threat1_ego == 1, 1, 0)


summary(lm(in_spokeref ~ both_5 + both_4  + both_3 + both_2 , data = L)) #it's matching on highest cat

m_views <- lm(in_spokeref ~ both_5 + both_4  + both_3 + both_2 + both_1, data = L)


## Appendix table A3:
stargazer(m_views)



## Let's repeat for importance

L$both_5_imp <- ifelse(L$interest_homoph == 1 & L$important_surr_refug_ego == 5, 1, 0)
L$both_4_imp <- ifelse(L$interest_homoph == 1 & L$important_surr_refug_ego == 4, 1, 0)
L$both_3_imp <- ifelse(L$interest_homoph == 1 & L$important_surr_refug_ego == 3, 1, 0)
L$both_2_imp <- ifelse(L$interest_homoph == 1 & L$important_surr_refug_ego == 2, 1, 0)
L$both_1_imp <- ifelse(L$interest_homoph == 1 & L$important_surr_refug_ego == 1, 1, 0)


## There are none that match on the lowest importance of 5 
## Once again, it's matching on most importance, not least importance, that matters.
summary(lm(in_spokeref ~  both_4_imp + both_3_imp + both_2_imp + both_1_imp , data = L)) #it's matching on highest cat

m_imp <- lm(in_spokeref ~  both_4_imp + both_3_imp + both_2_imp + both_1_imp , data = L)

stargazer(m_imp)


## We can add the highest categories for both to our QAP

L$highimp_homoph <- ifelse(L$interest_homoph == 1 & L$important_surr_refug_ego == 1, 1, 0)
L$posview_homoph <- ifelse(L$views_homoph == 1 & L$threat1_ego == 5, 1, 0)


## And the p-values:
summary(lm(in_spokeref ~ highimp_homoph, data = L))
summary(lm(in_spokeref ~ posview_homoph, data = L))

summary(lm(highimp_homoph ~ in_spokeref, data = L)) #.01
summary(lm(posview_homoph ~ in_spokeref, data = L)) #.00


## And these values, for addition to Table 6:

## let's remake these indicators for the L dataset

srsL <- which(L$in_spokeref == 1 & L$in_union == 1)

## Make an index to identify the rows in Links that are NOT in the spokeref network
## but ARE in the union social network
notsrsL <- which(L$in_spokeref == 0 & L$in_union == 1)

mean(L$highimp_homoph[srsL], na.rm=T) #.50
mean(L$highimp_homoph[notsrsL], na.rm=T) #.45

mean(L$posview_homoph[srsL], na.rm=T) #.24
mean(L$posview_homoph[notsrsL], na.rm=T) #.20

## And the p-values

t.test(L$highimp_homoph[srsL], L$highimp_homoph[notsrsL]) #.01
t.test(L$posview_homoph[srsL], L$posview_homoph[notsrsL]) #.00


############
## Quadratic Assignment Procedure for Table 7
## Exploring significance of coefficient alter attribute or homophily regressed
## on whether the link was used (difference in means between link used v. not 
## in terms of these attributes). 

## Once again, we permute attributes among the nodes within each village and 
## then re-run these regressions.  These once again hold the network structure
## and link use fixed and shuffle the attributes among a village's nodes.

nsims <- 1000
set.seed(20240419)

numvlgs <- length(unique(D$village_number))

sim_reg_coefs_3 <- matrix(data = NA, nrow = nsims, ncol = 12)
colnames(sim_reg_coefs_3) <- c("alter_ref", "alter_knows", "alter_farmer", "alter_views", "alter_int",
                               "relig_homoph", "lang_homoph", "ref_homoph", "views_homoph", "int_homoph",
                               "posview_homoph", "highimp_homoph")

L_sm<- L %>% select(ego_hh, alter_hh, village_number, in_spokeref)

## We also need to reconstruct the relevant node-level dataset with attributes

atts <- c("been_refugee",
          "refugee_meet",
          "threat1", 
          "farmer",
          "resp_language",
          "religion",
          "important_surr_refug")

## Now we need to create a dataset that includes the household number,
## village number, and these attributes.  

## First grab these attributes in their own dataset
egoatts <- D[,which(names(D) %in% atts)]

## Now add "_ego" to the names of these variables
names(egoatts) <- paste(names(egoatts), "ego", sep = "_")

## And add the household number which we call ego_hh
egoatts$ego_hh <- D$hhnum
## And the village number
egoatts$village_number <- D$village_number

## Repeat this same set of steps, this time calling the household number alter_hh
## adding "_alter" to the variable names, and merging by alter_hh and village number

alteratts <- D[,which(names(D) %in% atts)]

## Now add "_ego" to the names of these variables
names(alteratts) <- paste(names(alteratts), "alter", sep = "_")

## And add the household number which we'll call alter_hh
alteratts$alter_hh <- D$hhnum
## And the village number
alteratts$village_number <- D$village_number

## Now egoatts and alteratts are our two attribute datasets




## Confirm that the two attrbiute dataframes are in the same order before beginning
all.equal(egoatts$village_number, alteratts$village_number) #True
all.equal(egoatts$ego_hh, alteratts$alter_hh) #True

for (i in 1:nsims){
  
  ## randomize the row order for the egoatts (which will be the same for the alteratts)
  ## but do so within villages
  ## then randomly reorder those rows within villages
  ## then merge both into the links dataset
  ## then run all the regressions and store the coefficients
  
  ## Make the attribute matrices to shuffle
  sim_ego_atts <- egoatts
  sim_alter_atts <- alteratts
  
  ## For each village, randomly shuffle the rows of the egoatts and alteratts matrix
  
  for(j in 1:numvlgs){
    
    ## Identify the rows of the attribute matrix that are part of village j
    thevlg <- which(sim_ego_atts$village_number == j)
    
    ## Shuffle these rows by using an index randomly shuffled
    
    shuffled_vlg <- sample(1:length(thevlg), length(thevlg), replace = FALSE)
    
    ## Now replace each of that village's attributes with the same that have been shuffled
    
    ## Been refugee:
    sim_ego_atts$been_refugee_ego[thevlg] <- sim_ego_atts$been_refugee_ego[thevlg][shuffled_vlg]
    sim_alter_atts$been_refugee_alter[thevlg] <- sim_alter_atts$been_refugee_alter[thevlg][shuffled_vlg]
    
    
    ## Refugee meet:
    sim_ego_atts$refugee_meet_ego[thevlg] <- sim_ego_atts$refugee_meet_ego[thevlg][shuffled_vlg]
    sim_alter_atts$refugee_meet_alter[thevlg] <- sim_alter_atts$refugee_meet_alter[thevlg][shuffled_vlg]
    
    ## Resp language:
    sim_ego_atts$resp_language_ego[thevlg] <- sim_ego_atts$resp_language_ego[thevlg][shuffled_vlg]
    sim_alter_atts$resp_language_alter[thevlg] <- sim_alter_atts$resp_language_alter[thevlg][shuffled_vlg]
    
    ## Religion:
    sim_ego_atts$religion_ego[thevlg] <- sim_ego_atts$religion_ego[thevlg][shuffled_vlg]
    sim_alter_atts$religion_alter[thevlg] <-  sim_alter_atts$religion_alter[thevlg][shuffled_vlg]
    
    ## Farmer:
    sim_ego_atts$farmer_ego[thevlg] <- sim_ego_atts$farmer_ego[thevlg][shuffled_vlg]
    sim_alter_atts$farmer_alter[thevlg] <- sim_alter_atts$farmer_alter[thevlg][shuffled_vlg]
    
    ## Interest:
    sim_ego_atts$important_surr_refug_ego[thevlg] <- sim_ego_atts$important_surr_refug_ego[thevlg][shuffled_vlg]
    sim_alter_atts$important_surr_refug_alter[thevlg] <- sim_alter_atts$important_surr_refug_alter[thevlg][shuffled_vlg]
    
    ## Views: 
    sim_ego_atts$threat1_ego[thevlg] <- sim_ego_atts$threat1_ego[thevlg][shuffled_vlg]
    sim_alter_atts$threat1_alter[thevlg] <- sim_alter_atts$threat1_alter[thevlg][shuffled_vlg]
    
  }
  
  ## Now sim_ego_atts and sim_alter_atts has been shuffled by village on all 
  ## relevant attributes
  
  ## Now we need to merge these into our link dataset by first merging with L_sm
  
  L_sim <- merge(L_sm, sim_ego_atts, by = c("ego_hh", "village_number"), all.x = TRUE, all.y = FALSE)
  
  ## And then merging the output, L_sim, with the alter atts
  
  L_sim <- merge(L_sim, sim_alter_atts, by = c("alter_hh", "village_number"), all.x = TRUE, all.y = FALSE)
  
  ## Now we add our 7 homophily values to L_sim
  
  L_sim$relig_homoph <- ifelse(L_sim$religion_ego == L_sim$religion_alter, 1, 0)
  L_sim$lang_homoph <- ifelse(L_sim$resp_language_ego == L_sim$resp_language_alter, 1, 0)
  L_sim$ref_homoph <- ifelse(L_sim$been_refugee_ego == L_sim$been_refugee_alter, 1, 0)
  L_sim$views_homoph <- ifelse(L_sim$threat1_ego == L_sim$threat1_alter, 1, 0)
  L_sim$interest_homoph <- ifelse(L_sim$important_surr_refug_ego == L_sim$important_surr_refug_alter, 1, 0)
  L_sim$highimp_homoph <- ifelse(L_sim$interest_homoph == 1 & L_sim$important_surr_refug_ego == 1, 1, 0)
  L_sim$posview_homoph <- ifelse(L_sim$views_homoph == 1 & L_sim$threat1_ego == 5, 1, 0)
  
  
  ## Now run the twelve simple regressions and store the coefficients from each
  sim_reg_coefs_3[i,1] <- coef(lm(been_refugee_alter ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,2] <- coef(lm(refugee_meet_alter ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,3] <- coef(lm(farmer_alter ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,4] <- coef(lm(threat1_alter ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,5] <- coef(lm(important_surr_refug_alter ~ in_spokeref, data = L_sim))[2]

  ## And the homophily ones
  sim_reg_coefs_3[i,6] <- coef(lm(relig_homoph ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,7] <- coef(lm(lang_homoph ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,8] <- coef(lm(ref_homoph ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,9] <- coef(lm(views_homoph ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,10] <- coef(lm(interest_homoph ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,11] <- coef(lm(posview_homoph ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_3[i,12] <- coef(lm(highimp_homoph ~ in_spokeref, data = L_sim))[2]
  
  if(i %% 100 == 0) print(i)
  
}

## Next, store the correct reference coefficient



coef_beenref2 <- coef(lm(been_refugee_alter ~ in_spokeref, data = L))[2]
coef_metref2 <- coef(lm(refugee_meet_alter ~ in_spokeref, data = L))[2]
coef_farmer2 <- coef(lm(farmer_alter ~ in_spokeref, data = L))[2]
coef_views2 <- coef(lm(threat1_alter ~ in_spokeref, data = L))[2]
coef_int2 <- coef(lm(important_surr_refug_alter ~ in_spokeref, data = L))[2]

coef_relig_h2 <- coef(lm(relig_homoph ~ in_spokeref, data = L))[2]
coef_lang_h2 <- coef(lm(lang_homoph ~ in_spokeref, data = L))[2]
coef_beenref_h2 <- coef(lm(ref_homoph ~ in_spokeref, data = L))[2]
coef_views_h2 <- coef(lm(views_homoph ~ in_spokeref, data = L))[2]
coef_int_h2 <- coef(lm(interest_homoph ~ in_spokeref, data = L))[2]

coef_posviews_h2 <- coef(lm(posview_homoph ~ in_spokeref, data = L))[2]
coef_highint_h2 <- coef(lm(highimp_homoph ~ in_spokeref, data = L))[2]


## And now assess significance
sum(sim_reg_coefs_3[,1] >= coef_beenref2) #733
sum(sim_reg_coefs_3[,2] >= coef_metref2) #155
sum(sim_reg_coefs_3[,3] >= coef_farmer2) #360
sum(sim_reg_coefs_3[,4] >= coef_views2) #490
sum(sim_reg_coefs_3[,5] <= coef_int2) #1 THIS ONE IS A NEGATIVE COEFFICIENT

sum(sim_reg_coefs_3[,6] >= coef_relig_h2) #63
sum(sim_reg_coefs_3[,7] >= coef_lang_h2) #47
sum(sim_reg_coefs_3[,8] >= coef_beenref_h2) #303
sum(sim_reg_coefs_3[,9] >= coef_views_h2) #129
sum(sim_reg_coefs_3[,10] >= coef_int_h2) #30

sum(sim_reg_coefs_3[,11] >= coef_posviews_h2) #42
sum(sim_reg_coefs_3[,12] >= coef_highint_h2) #68

## And we can generate the two-sided QAP scores with:

sum(abs(sim_reg_coefs_3[,1]) >= abs(coef_beenref2))/1000 #.798
sum(abs(sim_reg_coefs_3[,2]) >= abs(coef_metref2))/1000 #.557
sum(abs(sim_reg_coefs_3[,3]) >= abs(coef_farmer2))/1000 #.672
sum(abs(sim_reg_coefs_3[,4]) >= abs(coef_views2))/1000 #.527
sum(abs(sim_reg_coefs_3[,5]) >= abs(coef_int2))/1000 #.001

sum(abs(sim_reg_coefs_3[,6]) >= abs(coef_relig_h2))/1000 #.20
sum(abs(sim_reg_coefs_3[,7]) >= abs(coef_lang_h2))/1000 #.747
sum(abs(sim_reg_coefs_3[,8]) >= abs(coef_beenref_h2))/1000 #.727
sum(abs(sim_reg_coefs_3[,9]) >= abs(coef_views_h2))/1000 #.149
sum(abs(sim_reg_coefs_3[,10]) >= abs(coef_int_h2))/1000 #.037

sum(abs(sim_reg_coefs_3[,11]) >= abs(coef_posviews_h2))/1000 #.042
sum(abs(sim_reg_coefs_3[,12]) >= abs(coef_highint_h2))/1000 #.078

## These go in Table 7

## Now let's plot some from this set for the paper
## We'll need the regression coefficient matrix to be a dataframe

sim_regs <- as.data.frame(sim_reg_coefs_3)

sim_regs %>%
  ggplot()+
  geom_density(aes(alter_ref), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_beenref2, linetype = "dashed")+
  labs(title = "QAP, Alter Was Refugee", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg1.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(alter_knows), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_metref2, linetype = "dashed")+
  labs(title = "QAP, Alter Knows Refugee", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg2.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(alter_farmer), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_farmer2, linetype = "dashed")+
  labs(title = "QAP, Alter Farmer", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg3.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(alter_views), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_views2, linetype = "dashed")+
  labs(title = "QAP, Alter's Views", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg4.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(alter_int), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_int2, linetype = "dashed")+
  labs(title = "QAP, Alter's Interest", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg5.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(relig_homoph), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_relig_h2, linetype = "dashed")+
  labs(title = "QAP, Religion Homophily", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg6.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(lang_homoph), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_lang_h2, linetype = "dashed")+
  labs(title = "QAP, Language Homophily", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg7.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(ref_homoph), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_beenref_h2, linetype = "dashed")+
  labs(title = "QAP, Refugee Status Homophily", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg8.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(views_homoph), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_views_h2, linetype = "dashed")+
  labs(title = "QAP, Refugee Views Homophily", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg9.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(int_homoph), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_int_h2, linetype = "dashed")+
  labs(title = "QAP, Interest Homophily", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg10.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(posview_homoph), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_posviews_h2, linetype = "dashed")+
  labs(title = "QAP, Positive Refugee Views Homophily", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg11.pdf

sim_regs %>%
  ggplot()+
  geom_density(aes(highimp_homoph), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_highint_h2, linetype = "dashed")+
  labs(title = "QAP, High Interest Homophily", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg12.pdf






#####
## Next, perform the same comparisons, this time looking only at the links used
## and contrasting between treated and control
## April 2024: add four rows for positive views homophily, high interest homophily,
## in multiple layers, and count of layers

## Let's add the four relevant variables to the Links dataset too (not just L as above)

Links$views_homoph <- ifelse(Links$threat1_ego == Links$threat1_alter, 1, 0)
Links$interest_homoph <- ifelse(Links$important_surr_refug_ego == Links$important_surr_refug_alter, 1, 0)
Links$highimp_homoph <- ifelse(Links$interest_homoph == 1 & Links$important_surr_refug_ego == 1, 1, 0)
Links$posview_homoph <- ifelse(Links$views_homoph == 1 & Links$threat1_ego == 5, 1, 0)
Links$multiplicity_count <- Links$in_meal + Links$in_visit + Links$in_rumor + Links$in_borrow
Links$multiplicity_indicator <- ifelse(Links$multiplicity_count > 1, 1, 0)


## Make an indicator for which links were used by a treated ego
Tsrs <- which(Links$in_spokeref == 1 & Links$trt_ego == 1)
## Make an indicator for which links were used by a control ego
Csrs <- which(Links$in_spokeref == 1 & Links$trt_ego == 0)

treateffect <- matrix(data = NA, nrow = 3, ncol = 15)
treateffect[1,1] <- length(Csrs) #number of control links used to discuss refs
treateffect[2,1] <- length(Tsrs) #number of treatment links used to discuss refs

## Are they seeking out someone with relevant experience?
#### Has been refugee?
treateffect[1,2] <- mean(Links$been_refugee_alter[Csrs], na.rm=TRUE) 
treateffect[2,2] <- mean(Links$been_refugee_alter[Tsrs], na.rm=TRUE)
treateffect[3,2] <- t.test(Links$been_refugee_alter[Csrs], Links$been_refugee_alter[Tsrs])$p.value


#### Knows refugee?
treateffect[1,3] <- mean(Links$refugee_meet_alter[Csrs], na.rm=T)
treateffect[2,3] <- mean(Links$refugee_meet_alter[Tsrs], na.rm=T)
treateffect[3,3] <- t.test(Links$refugee_meet_alter[Csrs], Links$refugee_meet_alter[Tsrs])$p.value


#### Farmer?
treateffect[1,4] <- mean(Links$farmer_alter[Csrs], na.rm=T)
treateffect[2,4] <- mean(Links$farmer_alter[Tsrs], na.rm=T)
treateffect[3,4] <- t.test(Links$farmer_alter[Csrs], Links$farmer_alter[Tsrs])$p.value

#### Feels a certain way about the topic?
treateffect[1,5] <- mean(Links$threat1_alter[Csrs], na.rm=T)
treateffect[2,5] <- mean(Links$threat1_alter[Tsrs], na.rm=T)
treateffect[3,5] <- t.test(Links$threat1_alter[Csrs], Links$threat1_alter[Tsrs])$p.value #not sig

#### Finds the topic to be important?

treateffect[1,6] <- mean(Links$important_surr_refug_alter[Csrs], na.rm=T)
treateffect[2,6] <- mean(Links$important_surr_refug_alter[Tsrs], na.rm=T)
treateffect[3,6] <- t.test(Links$important_surr_refug_alter[Csrs], Links$important_surr_refug_alter[Tsrs])$p.value #ok, importance of alter matters

## Is culturally similar to the respondent in terms of religion?

treateffect[1,7] <- mean(Links$religion_ego[Csrs] == Links$religion_alter[Csrs], na.rm=T)
treateffect[2,7] <- mean(Links$religion_ego[Tsrs] == Links$religion_alter[Tsrs], na.rm=T)
treateffect[3,7] <- t.test(Links$religion_ego[Csrs] == Links$religion_alter[Csrs],
                           Links$religion_ego[Tsrs] == Links$religion_alter[Tsrs])$p.value #not sig, close-ish

#### In terms of language?

treateffect[1,8] <- mean(Links$resp_language_ego[Csrs] == Links$resp_language_alter[Csrs], na.rm=T)
treateffect[2,8] <- mean(Links$resp_language_ego[Tsrs] == Links$resp_language_alter[Tsrs], na.rm=T)
treateffect[3,8] <- t.test(Links$resp_language_ego[Csrs] == Links$resp_language_alter[Csrs],
                           Links$resp_language_ego[Tsrs] == Links$resp_language_alter[Tsrs])$p.value

#### In terms of shared refugee experience?

treateffect[1,9] <- mean(Links$been_refugee_ego[Csrs] == Links$been_refugee_alter[Csrs], na.rm=T)
treateffect[2,9] <- mean(Links$been_refugee_ego[Tsrs] == Links$been_refugee_alter[Tsrs], na.rm=T)
treateffect[3,9] <- t.test(Links$been_refugee_ego[Csrs] == Links$been_refugee_alter[Csrs],
                           Links$been_refugee_ego[Tsrs] == Links$been_refugee_alter[Tsrs])$p.value #not homoph w.r.t. ref status

#### In terms of shared views on the topic?
treateffect[1,10] <- mean(Links$threat1_ego[Csrs] == Links$threat1_alter[Csrs], na.rm=T)
treateffect[2,10] <- mean(Links$threat1_ego[Tsrs] == Links$threat1_alter[Tsrs], na.rm=T)
treateffect[3,10] <- t.test(Links$threat1_ego[Csrs] == Links$threat1_alter[Csrs], Links$threat1_ego[Tsrs] == Links$threat1_alter[Tsrs])$p.value #sig at .1, seek out similarit

#### In terms of shared perception of importance?
treateffect[1,11] <- mean(Links$important_surr_refug_ego[Csrs] == Links$important_surr_refug_alter[Csrs], na.rm=T)
treateffect[2,11] <- mean(Links$important_surr_refug_ego[Tsrs] == Links$important_surr_refug_alter[Tsrs], na.rm=T)
treateffect[3,11] <- t.test(Links$important_surr_refug_ego[Csrs] == Links$important_surr_refug_alter[Csrs],
                            Links$important_surr_refug_ego[Tsrs] == Links$important_surr_refug_alter[Tsrs]  )$p.value #sig at .5


#### In terms of shared positive view on topic?
treateffect[1,12] <- mean(Links$posview_homoph[Csrs], na.rm=T)
treateffect[2,12] <- mean(Links$posview_homoph[Tsrs], na.rm=T)
treateffect[3,12] <- t.test(Links$posview_homoph[Csrs], Links$posview_homoph[Tsrs])$p.value

#### In terms of shared high interest in topic?
treateffect[1,13] <- mean(Links$highimp_homoph[Csrs], na.rm=T)
treateffect[2,13] <- mean(Links$highimp_homoph[Tsrs], na.rm=T)
treateffect[3,13] <- t.test(Links$highimp_homoph[Csrs], Links$highimp_homoph[Tsrs])$p.value

#### In terms of mutliplexity?
treateffect[1,14] <- mean(Links$multiplicity_indicator[Csrs], na.rm=T)
treateffect[2,14] <- mean(Links$multiplicity_indicator[Tsrs], na.rm=T)
treateffect[3,14] <- t.test(Links$multiplicity_indicator[Csrs], Links$multiplicity_indicator[Tsrs])$p.value


#### In terms of a count of layers?
treateffect[1,15] <- mean(Links$multiplicity_count[Csrs], na.rm=T)
treateffect[2,15] <- mean(Links$multiplicity_count[Tsrs], na.rm=T)
treateffect[3,15] <- t.test(Links$multiplicity_count[Csrs], Links$multiplicity_count[Tsrs])$p.value

colnames(treateffect) <- c("Links", "Alter was ref", "Alter knows ref", "Alter farmer", "Alter's views", "Alter's interest", "Relig homoph", "Language homoph", 
                           "Refugee status homoph", "Refugee views homoph", "Interest homoph", "Positive views homoph", "High interest homoph",
                           "In Multiple Layers", "Count of Layers")

rownames(treateffect) <- c("Control", "Treatment", "p-value")

## And build Paper Table 7
xtable(t(treateffect), digits = 2)



##########
## Now repeat, separating by network layer type
## to build Appendix Tables A3-A6

## --------------- Comparing ATE across Meal Links ---------------------------
Links_meal <- subset(Links, Links$in_meal == 1)
Links_visit <- subset(Links, Links$in_visit == 1)
Links_rumor <- subset(Links, Links$in_rumor == 1)
Links_borrow <- subset(Links, Links$in_borrow == 1)

## Meal. 
Tsrs_meal <- which(Links_meal$in_spokeref == 1 & Links_meal$trt_ego == 1)

Csrs_meal <- which(Links_meal$in_spokeref == 1 & Links_meal$trt_ego == 0)


treateffect_meal <- matrix(data = NA, nrow = 3, ncol = 15)
treateffect_meal[1,1] <- length(Csrs_meal) #number of control links used to discuss refs
treateffect_meal[2,1] <- length(Tsrs_meal) #number of treatment links used to discuss refs

## Are they seeking out someone with relevant experience?
#### Has been refugee?
treateffect_meal[1,2] <- mean(Links_meal$been_refugee_alter[Csrs_meal], na.rm=TRUE) 
treateffect_meal[2,2] <- mean(Links_meal$been_refugee_alter[Tsrs_meal], na.rm=TRUE)
treateffect_meal[3,2] <- t.test(Links_meal$been_refugee_alter[Csrs_meal], Links_meal$been_refugee_alter[Tsrs_meal])$p.value


#### Knows refugee?
treateffect_meal[1,3] <- mean(Links_meal$refugee_meet_alter[Csrs_meal], na.rm=T)
treateffect_meal[2,3] <- mean(Links_meal$refugee_meet_alter[Tsrs_meal], na.rm=T)
treateffect_meal[3,3] <- t.test(Links_meal$refugee_meet_alter[Csrs_meal], Links_meal$refugee_meet_alter[Tsrs_meal])$p.value


#### Farmer?
treateffect_meal[1,4] <- mean(Links_meal$farmer_alter[Csrs_meal], na.rm=T)
treateffect_meal[2,4] <- mean(Links_meal$farmer_alter[Tsrs_meal], na.rm=T)
treateffect_meal[3,4] <- t.test(Links_meal$farmer_alter[Csrs_meal], Links_meal$farmer_alter[Tsrs_meal])$p.value

#### Feels a certain way about the topic?
treateffect_meal[1,5] <- mean(Links_meal$threat1_alter[Csrs_meal], na.rm=T)
treateffect_meal[2,5] <- mean(Links_meal$threat1_alter[Tsrs_meal], na.rm=T)
treateffect_meal[3,5] <- t.test(Links_meal$threat1_alter[Csrs_meal], Links_meal$threat1_alter[Tsrs_meal])$p.value #not sig

#### Finds the topic to be important?

treateffect_meal[1,6] <- mean(Links_meal$important_surr_refug_alter[Csrs_meal], na.rm=T)
treateffect_meal[2,6] <- mean(Links_meal$important_surr_refug_alter[Tsrs_meal], na.rm=T)
treateffect_meal[3,6] <- t.test(Links_meal$important_surr_refug_alter[Csrs_meal], Links_meal$important_surr_refug_alter[Tsrs_meal])$p.value #ok, importance of alter matters

## Is culturally similar to the respondent in terms of religion?

treateffect_meal[1,7] <- mean(Links_meal$religion_ego[Csrs_meal] == Links_meal$religion_alter[Csrs_meal], na.rm=T)
treateffect_meal[2,7] <- mean(Links_meal$religion_ego[Tsrs_meal] == Links_meal$religion_alter[Tsrs_meal], na.rm=T)
treateffect_meal[3,7] <- t.test(Links_meal$religion_ego[Csrs_meal] == Links_meal$religion_alter[Csrs_meal],
                                Links_meal$religion_ego[Tsrs_meal] == Links_meal$religion_alter[Tsrs_meal])$p.value #not sig, close-ish

#### In terms of language?

treateffect_meal[1,8] <- mean(Links_meal$resp_language_ego[Csrs_meal] == Links_meal$resp_language_alter[Csrs_meal], na.rm=T)
treateffect_meal[2,8] <- mean(Links_meal$resp_language_ego[Tsrs_meal] == Links_meal$resp_language_alter[Tsrs_meal], na.rm=T)
treateffect_meal[3,8] <- t.test(Links_meal$resp_language_ego[Csrs_meal] == Links_meal$resp_language_alter[Csrs_meal],
                                Links_meal$resp_language_ego[Tsrs_meal] == Links_meal$resp_language_alter[Tsrs_meal])$p.value

#### In terms of shared refugee experience?

treateffect_meal[1,9] <- mean(Links_meal$been_refugee_ego[Csrs_meal] == Links_meal$been_refugee_alter[Csrs_meal], na.rm=T)
treateffect_meal[2,9] <- mean(Links_meal$been_refugee_ego[Tsrs_meal] == Links_meal$been_refugee_alter[Tsrs_meal], na.rm=T)
treateffect_meal[3,9] <- t.test(Links_meal$been_refugee_ego[Csrs_meal] == Links_meal$been_refugee_alter[Csrs_meal],
                                Links_meal$been_refugee_ego[Tsrs_meal] == Links_meal$been_refugee_alter[Tsrs_meal])$p.value #not homoph w.r.t. ref status

#### In terms of shared views on the topic?
treateffect_meal[1,10] <- mean(Links_meal$threat1_ego[Csrs_meal] == Links_meal$threat1_alter[Csrs_meal], na.rm=T)
treateffect_meal[2,10] <- mean(Links_meal$threat1_ego[Tsrs_meal] == Links_meal$threat1_alter[Tsrs_meal], na.rm=T)
treateffect_meal[3,10] <- t.test(Links_meal$threat1_ego[Csrs_meal] == Links_meal$threat1_alter[Csrs_meal], Links_meal$threat1_ego[Tsrs_meal] == Links_meal$threat1_alter[Tsrs_meal])$p.value #sig at .1, seek out similarit

#### In terms of shared perception of importance?
treateffect_meal[1,11] <- mean(Links_meal$important_surr_refug_ego[Csrs_meal] == Links_meal$important_surr_refug_alter[Csrs_meal], na.rm=T)
treateffect_meal[2,11] <- mean(Links_meal$important_surr_refug_ego[Tsrs_meal] == Links_meal$important_surr_refug_alter[Tsrs_meal], na.rm=T)
treateffect_meal[3,11] <- t.test(Links_meal$important_surr_refug_ego[Csrs_meal] == Links_meal$important_surr_refug_alter[Csrs_meal],
                                 Links_meal$important_surr_refug_ego[Tsrs_meal] == Links_meal$important_surr_refug_alter[Tsrs_meal]  )$p.value #sig at .5

#### In terms of shared positive view on topic?
treateffect_meal[1,12] <- mean(Links_meal$posview_homoph[Csrs_meal], na.rm=T)
treateffect_meal[2,12] <- mean(Links_meal$posview_homoph[Tsrs_meal], na.rm=T)
treateffect_meal[3,12] <- t.test(Links_meal$posview_homoph[Csrs_meal], Links_meal$posview_homoph[Tsrs_meal])$p.value

#### In terms of shared high interest in topic?
treateffect_meal[1,13] <- mean(Links_meal$highimp_homoph[Csrs_meal], na.rm=T)
treateffect_meal[2,13] <- mean(Links_meal$highimp_homoph[Tsrs_meal], na.rm=T)
treateffect_meal[3,13] <- t.test(Links_meal$highimp_homoph[Csrs_meal], Links_meal$highimp_homoph[Tsrs_meal])$p.value

#### In terms of mutliplexity?
treateffect_meal[1,14] <- mean(Links_meal$multiplicity_indicator[Csrs_meal], na.rm=T)
treateffect_meal[2,14] <- mean(Links_meal$multiplicity_indicator[Tsrs_meal], na.rm=T)
treateffect_meal[3,14] <- t.test(Links_meal$multiplicity_indicator[Csrs_meal], Links_meal$multiplicity_indicator[Tsrs_meal])$p.value


#### In terms of a count of layers?
treateffect_meal[1,15] <- mean(Links$multiplicity_count[Csrs], na.rm=T)
treateffect_meal[2,15] <- mean(Links$multiplicity_count[Tsrs], na.rm=T)
treateffect_meal[3,15] <- t.test(Links$multiplicity_count[Csrs], Links$multiplicity_count[Tsrs])$p.value

colnames(treateffect_meal) <- c("Links", "Alter was ref", "Alter knows ref", "Alter farmer", "Alter's views", "Alter's interest", "Relig homoph", "Language homoph", 
                           "Refugee status homoph", "Refugee views homoph", "Interest homoph", "Positive views homoph", "High interest homoph",
                           "In Multiple Layers", "Count of Layers")

rownames(treateffect_meal) <- c("Control", "Treatment", "p-value")



## SI Table A3
xtable(t(treateffect_meal), digits = 2)

## --------------- Comparing ATE across Visit Links ---------------------------
## Visit.  
Tsrs_visit <- which(Links_visit$in_spokeref == 1 & Links_visit$trt_ego == 1)

Csrs_visit <- which(Links_visit$in_spokeref == 1 & Links_visit$trt_ego == 0)


treateffect_visit <- matrix(data = NA, nrow = 3, ncol = 15)
treateffect_visit[1,1] <- length(Csrs_visit) #number of control links used to discuss refs
treateffect_visit[2,1] <- length(Tsrs_visit) #number of treatment links used to discuss refs

## Are they seeking out someone with relevant experience?
#### Has been refugee?
treateffect_visit[1,2] <- mean(Links_visit$been_refugee_alter[Csrs_visit], na.rm=TRUE) 
treateffect_visit[2,2] <- mean(Links_visit$been_refugee_alter[Tsrs_visit], na.rm=TRUE)
treateffect_visit[3,2] <- t.test(Links_visit$been_refugee_alter[Csrs_visit], Links_visit$been_refugee_alter[Tsrs_visit])$p.value


#### Knows refugee?
treateffect_visit[1,3] <- mean(Links_visit$refugee_meet_alter[Csrs_visit], na.rm=T)
treateffect_visit[2,3] <- mean(Links_visit$refugee_meet_alter[Tsrs_visit], na.rm=T)
treateffect_visit[3,3] <- t.test(Links_visit$refugee_meet_alter[Csrs_visit], Links_visit$refugee_meet_alter[Tsrs_visit])$p.value


#### Farmer?
treateffect_visit[1,4] <- mean(Links_visit$farmer_alter[Csrs_visit], na.rm=T)
treateffect_visit[2,4] <- mean(Links_visit$farmer_alter[Tsrs_visit], na.rm=T)
treateffect_visit[3,4] <- t.test(Links_visit$farmer_alter[Csrs_visit], Links_visit$farmer_alter[Tsrs_visit])$p.value

#### Feels a certain way about the topic?
treateffect_visit[1,5] <- mean(Links_visit$threat1_alter[Csrs_visit], na.rm=T)
treateffect_visit[2,5] <- mean(Links_visit$threat1_alter[Tsrs_visit], na.rm=T)
treateffect_visit[3,5] <- t.test(Links_visit$threat1_alter[Csrs_visit], Links_visit$threat1_alter[Tsrs_visit])$p.value #not sig

#### Finds the topic to be important?

treateffect_visit[1,6] <- mean(Links_visit$important_surr_refug_alter[Csrs_visit], na.rm=T)
treateffect_visit[2,6] <- mean(Links_visit$important_surr_refug_alter[Tsrs_visit], na.rm=T)
treateffect_visit[3,6] <- t.test(Links_visit$important_surr_refug_alter[Csrs_visit], Links_visit$important_surr_refug_alter[Tsrs_visit])$p.value #ok, importance of alter matters

## Is culturally similar to the respondent in terms of religion?

treateffect_visit[1,7] <- mean(Links_visit$religion_ego[Csrs_visit] == Links_visit$religion_alter[Csrs_visit], na.rm=T)
treateffect_visit[2,7] <- mean(Links_visit$religion_ego[Tsrs_visit] == Links_visit$religion_alter[Tsrs_visit], na.rm=T)
treateffect_visit[3,7] <- t.test(Links_visit$religion_ego[Csrs_visit] == Links_visit$religion_alter[Csrs_visit],
                                 Links_visit$religion_ego[Tsrs_visit] == Links_visit$religion_alter[Tsrs_visit])$p.value #not sig, close-ish

#### In terms of language?

treateffect_visit[1,8] <- mean(Links_visit$resp_language_ego[Csrs_visit] == Links_visit$resp_language_alter[Csrs_visit], na.rm=T)
treateffect_visit[2,8] <- mean(Links_visit$resp_language_ego[Tsrs_visit] == Links_visit$resp_language_alter[Tsrs_visit], na.rm=T)
treateffect_visit[3,8] <- t.test(Links_visit$resp_language_ego[Csrs_visit] == Links_visit$resp_language_alter[Csrs_visit],
                                 Links_visit$resp_language_ego[Tsrs_visit] == Links_visit$resp_language_alter[Tsrs_visit])$p.value

#### In terms of shared refugee experience?

treateffect_visit[1,9] <- mean(Links_visit$been_refugee_ego[Csrs_visit] == Links_visit$been_refugee_alter[Csrs_visit], na.rm=T)
treateffect_visit[2,9] <- mean(Links_visit$been_refugee_ego[Tsrs_visit] == Links_visit$been_refugee_alter[Tsrs_visit], na.rm=T)
treateffect_visit[3,9] <- t.test(Links_visit$been_refugee_ego[Csrs_visit] == Links_visit$been_refugee_alter[Csrs_visit],
                                 Links_visit$been_refugee_ego[Tsrs_visit] == Links_visit$been_refugee_alter[Tsrs_visit])$p.value #not homoph w.r.t. ref status

#### In terms of shared views on the topic?
treateffect_visit[1,10] <- mean(Links_visit$threat1_ego[Csrs_visit] == Links_visit$threat1_alter[Csrs_visit], na.rm=T)
treateffect_visit[2,10] <- mean(Links_visit$threat1_ego[Tsrs_visit] == Links_visit$threat1_alter[Tsrs_visit], na.rm=T)
treateffect_visit[3,10] <- t.test(Links_visit$threat1_ego[Csrs_visit] == Links_visit$threat1_alter[Csrs_visit], 
                                  Links_visit$threat1_ego[Tsrs_visit] == Links_visit$threat1_alter[Tsrs_visit])$p.value #sig at .1, seek out similarit

#### In terms of shared perception of importance?
treateffect_visit[1,11] <- mean(Links_visit$important_surr_refug_ego[Csrs_visit] == Links_visit$important_surr_refug_alter[Csrs_visit], na.rm=T)
treateffect_visit[2,11] <- mean(Links_visit$important_surr_refug_ego[Tsrs_visit] == Links_visit$important_surr_refug_alter[Tsrs_visit], na.rm=T)
treateffect_visit[3,11] <- t.test(Links_visit$important_surr_refug_ego[Csrs_visit] == Links_visit$important_surr_refug_alter[Csrs_visit],
                                  Links_visit$important_surr_refug_ego[Tsrs_visit] == Links_visit$important_surr_refug_alter[Tsrs_visit]  )$p.value #sig at .5

#### In terms of shared positive view on topic?
treateffect_visit[1,12] <- mean(Links_visit$posview_homoph[Csrs_visit], na.rm=T)
treateffect_visit[2,12] <- mean(Links_visit$posview_homoph[Tsrs_visit], na.rm=T)
treateffect_visit[3,12] <- t.test(Links_visit$posview_homoph[Csrs_visit], Links_visit$posview_homoph[Tsrs_visit])$p.value

#### In terms of shared high interest in topic?
treateffect_visit[1,13] <- mean(Links_visit$highimp_homoph[Csrs_visit], na.rm=T)
treateffect_visit[2,13] <- mean(Links_visit$highimp_homoph[Tsrs_visit], na.rm=T)
treateffect_visit[3,13] <- t.test(Links_visit$highimp_homoph[Csrs_visit], Links_visit$highimp_homoph[Tsrs_visit])$p.value

#### In terms of mutliplexity?
treateffect_visit[1,14] <- mean(Links_visit$multiplicity_indicator[Csrs_visit], na.rm=T)
treateffect_visit[2,14] <- mean(Links_visit$multiplicity_indicator[Tsrs_visit], na.rm=T)
treateffect_visit[3,14] <- t.test(Links_visit$multiplicity_indicator[Csrs_visit], Links_visit$multiplicity_indicator[Tsrs_visit])$p.value


#### In terms of a count of layers?
treateffect_visit[1,15] <- mean(Links_visit$multiplicity_count[Csrs_visit], na.rm=T)
treateffect_visit[2,15] <- mean(Links_visit$multiplicity_count[Tsrs_visit], na.rm=T)
treateffect_visit[3,15] <- t.test(Links_visit$multiplicity_count[Csrs_visit], Links_visit$multiplicity_count[Tsrs_visit])$p.value

colnames(treateffect_visit) <- c("Links", "Alter was ref", "Alter knows ref", "Alter farmer", "Alter's views", "Alter's interest", "Relig homoph", "Language homoph", 
                                "Refugee status homoph", "Refugee views homoph", "Interest homoph", "Positive views homoph", "High interest homoph",
                                "In Multiple Layers", "Count of Layers")

rownames(treateffect_visit) <- c("Control", "Treatment", "p-value")




## SI Table A4
xtable(t(treateffect_visit), digits = 2)

## --------------- Comparing ATE across Rumor Links ---------------------------
## Rumor  
Tsrs_rumor <- which(Links_rumor$in_spokeref == 1 & Links_rumor$trt_ego == 1)

Csrs_rumor <- which(Links_rumor$in_spokeref == 1 & Links_rumor$trt_ego == 0)


treateffect_rumor <- matrix(data = NA, nrow = 3, ncol = 15)
treateffect_rumor[1,1] <- length(Csrs_rumor) #number of control links used to discuss refs
treateffect_rumor[2,1] <- length(Tsrs_rumor) #number of treatment links used to discuss refs

## Are they seeking out someone with relevant experience?
#### Has been refugee?
treateffect_rumor[1,2] <- mean(Links_rumor$been_refugee_alter[Csrs_rumor], na.rm=TRUE) 
treateffect_rumor[2,2] <- mean(Links_rumor$been_refugee_alter[Tsrs_rumor], na.rm=TRUE)
treateffect_rumor[3,2] <- t.test(Links_rumor$been_refugee_alter[Csrs_rumor], Links_rumor$been_refugee_alter[Tsrs_rumor])$p.value


#### Knows refugee?
treateffect_rumor[1,3] <- mean(Links_rumor$refugee_meet_alter[Csrs_rumor], na.rm=T)
treateffect_rumor[2,3] <- mean(Links_rumor$refugee_meet_alter[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,3] <- t.test(Links_rumor$refugee_meet_alter[Csrs_rumor], Links_rumor$refugee_meet_alter[Tsrs_rumor])$p.value


#### Farmer?
treateffect_rumor[1,4] <- mean(Links_rumor$farmer_alter[Csrs_rumor], na.rm=T)
treateffect_rumor[2,4] <- mean(Links_rumor$farmer_alter[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,4] <- t.test(Links_rumor$farmer_alter[Csrs_rumor], Links_rumor$farmer_alter[Tsrs_rumor])$p.value

#### Feels a certain way about the topic?
treateffect_rumor[1,5] <- mean(Links_rumor$threat1_alter[Csrs_rumor], na.rm=T)
treateffect_rumor[2,5] <- mean(Links_rumor$threat1_alter[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,5] <- t.test(Links_rumor$threat1_alter[Csrs_rumor], Links_rumor$threat1_alter[Tsrs_rumor])$p.value #not sig

#### Finds the topic to be important?

treateffect_rumor[1,6] <- mean(Links_rumor$important_surr_refug_alter[Csrs_rumor], na.rm=T)
treateffect_rumor[2,6] <- mean(Links_rumor$important_surr_refug_alter[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,6] <- t.test(Links_rumor$important_surr_refug_alter[Csrs_rumor], Links_rumor$important_surr_refug_alter[Tsrs_rumor])$p.value #ok, importance of alter matters

## Is culturally similar to the respondent in terms of religion?

treateffect_rumor[1,7] <- mean(Links_rumor$religion_ego[Csrs_rumor] == Links_rumor$religion_alter[Csrs_rumor], na.rm=T)
treateffect_rumor[2,7] <- mean(Links_rumor$religion_ego[Tsrs_rumor] == Links_rumor$religion_alter[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,7] <- t.test(Links_rumor$religion_ego[Csrs_rumor] == Links_rumor$religion_alter[Csrs_rumor],
                                 Links_rumor$religion_ego[Tsrs_rumor] == Links_rumor$religion_alter[Tsrs_rumor])$p.value #not sig, close-ish

#### In terms of language?

treateffect_rumor[1,8] <- mean(Links_rumor$resp_language_ego[Csrs_rumor] == Links_rumor$resp_language_alter[Csrs_rumor], na.rm=T)
treateffect_rumor[2,8] <- mean(Links_rumor$resp_language_ego[Tsrs_rumor] == Links_rumor$resp_language_alter[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,8] <- t.test(Links_rumor$resp_language_ego[Csrs_rumor] == Links_rumor$resp_language_alter[Csrs_rumor],
                                 Links_rumor$resp_language_ego[Tsrs_rumor] == Links_rumor$resp_language_alter[Tsrs_rumor])$p.value

#### In terms of shared refugee experience?

treateffect_rumor[1,9] <- mean(Links_rumor$been_refugee_ego[Csrs_rumor] == Links_rumor$been_refugee_alter[Csrs_rumor], na.rm=T)
treateffect_rumor[2,9] <- mean(Links_rumor$been_refugee_ego[Tsrs_rumor] == Links_rumor$been_refugee_alter[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,9] <- t.test(Links_rumor$been_refugee_ego[Csrs_rumor] == Links_rumor$been_refugee_alter[Csrs_rumor],
                                 Links_rumor$been_refugee_ego[Tsrs_rumor] == Links_rumor$been_refugee_alter[Tsrs_rumor])$p.value #not homoph w.r.t. ref status

#### In terms of shared views on the topic?
treateffect_rumor[1,10] <- mean(Links_rumor$threat1_ego[Csrs_rumor] == Links_rumor$threat1_alter[Csrs_rumor], na.rm=T)
treateffect_rumor[2,10] <- mean(Links_rumor$threat1_ego[Tsrs_rumor] == Links_rumor$threat1_alter[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,10] <- t.test(Links_rumor$threat1_ego[Csrs_rumor] == Links_rumor$threat1_alter[Csrs_rumor], 
                                  Links_rumor$threat1_ego[Tsrs_rumor] == Links_rumor$threat1_alter[Tsrs_rumor])$p.value #sig at .1, seek out similarit

#### In terms of shared perception of importance?
treateffect_rumor[1,11] <- mean(Links_rumor$important_surr_refug_ego[Csrs_rumor] == Links_rumor$important_surr_refug_alter[Csrs_rumor], na.rm=T)
treateffect_rumor[2,11] <- mean(Links_rumor$important_surr_refug_ego[Tsrs_rumor] == Links_rumor$important_surr_refug_alter[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,11] <- t.test(Links_rumor$important_surr_refug_ego[Csrs_rumor] == Links_rumor$important_surr_refug_alter[Csrs_rumor],
                                  Links_rumor$important_surr_refug_ego[Tsrs_rumor] == Links_rumor$important_surr_refug_alter[Tsrs_rumor]  )$p.value #sig at .5


#### In terms of shared positive view on topic?
treateffect_rumor[1,12] <- mean(Links_rumor$posview_homoph[Csrs_rumor], na.rm=T)
treateffect_rumor[2,12] <- mean(Links_rumor$posview_homoph[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,12] <- t.test(Links_rumor$posview_homoph[Csrs_rumor], Links_rumor$posview_homoph[Tsrs_rumor])$p.value

#### In terms of shared high interest in topic?
treateffect_rumor[1,13] <- mean(Links_rumor$highimp_homoph[Csrs_rumor], na.rm=T)
treateffect_rumor[2,13] <- mean(Links_rumor$highimp_homoph[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,13] <- t.test(Links_rumor$highimp_homoph[Csrs_rumor], Links_rumor$highimp_homoph[Tsrs_rumor])$p.value

#### In terms of mutliplexity?
treateffect_rumor[1,14] <- mean(Links_rumor$multiplicity_indicator[Csrs_rumor], na.rm=T)
treateffect_rumor[2,14] <- mean(Links_rumor$multiplicity_indicator[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,14] <- t.test(Links_rumor$multiplicity_indicator[Csrs_rumor], Links_rumor$multiplicity_indicator[Tsrs_rumor])$p.value


#### In terms of a count of layers?
treateffect_rumor[1,15] <- mean(Links_rumor$multiplicity_count[Csrs_rumor], na.rm=T)
treateffect_rumor[2,15] <- mean(Links_rumor$multiplicity_count[Tsrs_rumor], na.rm=T)
treateffect_rumor[3,15] <- t.test(Links_rumor$multiplicity_count[Csrs_rumor], Links_rumor$multiplicity_count[Tsrs_rumor])$p.value

colnames(treateffect_rumor) <- c("Links", "Alter was ref", "Alter knows ref", "Alter farmer", "Alter's views", "Alter's interest", "Relig homoph", "Language homoph", 
                                 "Refugee status homoph", "Refugee views homoph", "Interest homoph", "Positive views homoph", "High interest homoph",
                                 "In Multiple Layers", "Count of Layers")

rownames(treateffect_rumor) <- c("Control", "Treatment", "p-value")


## And build Appendix Table A5
xtable(t(treateffect_rumor), digits = 2)

## --------------- Comparing ATE across Borrow Links ---------------------------
## Borrow  
Tsrs_borrow <- which(Links_borrow$in_spokeref == 1 & Links_borrow$trt_ego == 1)

Csrs_borrow <- which(Links_borrow$in_spokeref == 1 & Links_borrow$trt_ego == 0)


treateffect_borrow <- matrix(data = NA, nrow = 3, ncol = 15)
treateffect_borrow[1,1] <- length(Csrs_borrow) #number of control links used to discuss refs
treateffect_borrow[2,1] <- length(Tsrs_borrow) #number of treatment links used to discuss refs

## Are they seeking out someone with relevant experience?
#### Has been refugee?
treateffect_borrow[1,2] <- mean(Links_borrow$been_refugee_alter[Csrs_borrow], na.rm=TRUE) 
treateffect_borrow[2,2] <- mean(Links_borrow$been_refugee_alter[Tsrs_borrow], na.rm=TRUE)
treateffect_borrow[3,2] <- t.test(Links_borrow$been_refugee_alter[Csrs_borrow], Links_borrow$been_refugee_alter[Tsrs_borrow])$p.value


#### Knows refugee?
treateffect_borrow[1,3] <- mean(Links_borrow$refugee_meet_alter[Csrs_borrow], na.rm=T)
treateffect_borrow[2,3] <- mean(Links_borrow$refugee_meet_alter[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,3] <- t.test(Links_borrow$refugee_meet_alter[Csrs_borrow], Links_borrow$refugee_meet_alter[Tsrs_borrow])$p.value


#### Farmer?
treateffect_borrow[1,4] <- mean(Links_borrow$farmer_alter[Csrs_borrow], na.rm=T)
treateffect_borrow[2,4] <- mean(Links_borrow$farmer_alter[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,4] <- t.test(Links_borrow$farmer_alter[Csrs_borrow], Links_borrow$farmer_alter[Tsrs_borrow])$p.value

#### Feels a certain way about the topic?
treateffect_borrow[1,5] <- mean(Links_borrow$threat1_alter[Csrs_borrow], na.rm=T)
treateffect_borrow[2,5] <- mean(Links_borrow$threat1_alter[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,5] <- t.test(Links_borrow$threat1_alter[Csrs_borrow], Links_borrow$threat1_alter[Tsrs_borrow])$p.value #not sig

#### Finds the topic to be important?

treateffect_borrow[1,6] <- mean(Links_borrow$important_surr_refug_alter[Csrs_borrow], na.rm=T)
treateffect_borrow[2,6] <- mean(Links_borrow$important_surr_refug_alter[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,6] <- t.test(Links_borrow$important_surr_refug_alter[Csrs_borrow], Links_borrow$important_surr_refug_alter[Tsrs_borrow])$p.value #ok, importance of alter matters

## Is culturally similar to the respondent in terms of religion?

treateffect_borrow[1,7] <- mean(Links_borrow$religion_ego[Csrs_borrow] == Links_borrow$religion_alter[Csrs_borrow], na.rm=T)
treateffect_borrow[2,7] <- mean(Links_borrow$religion_ego[Tsrs_borrow] == Links_borrow$religion_alter[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,7] <- t.test(Links_borrow$religion_ego[Csrs_borrow] == Links_borrow$religion_alter[Csrs_borrow],
                                  Links_borrow$religion_ego[Tsrs_borrow] == Links_borrow$religion_alter[Tsrs_borrow])$p.value #not sig, close-ish

#### In terms of language?

treateffect_borrow[1,8] <- mean(Links_borrow$resp_language_ego[Csrs_borrow] == Links_borrow$resp_language_alter[Csrs_borrow], na.rm=T)
treateffect_borrow[2,8] <- mean(Links_borrow$resp_language_ego[Tsrs_borrow] == Links_borrow$resp_language_alter[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,8] <- t.test(Links_borrow$resp_language_ego[Csrs_borrow] == Links_borrow$resp_language_alter[Csrs_borrow],
                                  Links_borrow$resp_language_ego[Tsrs_borrow] == Links_borrow$resp_language_alter[Tsrs_borrow])$p.value

#### In terms of shared refugee experience?

treateffect_borrow[1,9] <- mean(Links_borrow$been_refugee_ego[Csrs_borrow] == Links_borrow$been_refugee_alter[Csrs_borrow], na.rm=T)
treateffect_borrow[2,9] <- mean(Links_borrow$been_refugee_ego[Tsrs_borrow] == Links_borrow$been_refugee_alter[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,9] <- t.test(Links_borrow$been_refugee_ego[Csrs_borrow] == Links_borrow$been_refugee_alter[Csrs_borrow],
                                  Links_borrow$been_refugee_ego[Tsrs_borrow] == Links_borrow$been_refugee_alter[Tsrs_borrow])$p.value #not homoph w.r.t. ref status

#### In terms of shared views on the topic?
treateffect_borrow[1,10] <- mean(Links_borrow$threat1_ego[Csrs_borrow] == Links_borrow$threat1_alter[Csrs_borrow], na.rm=T)
treateffect_borrow[2,10] <- mean(Links_borrow$threat1_ego[Tsrs_borrow] == Links_borrow$threat1_alter[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,10] <- t.test(Links_borrow$threat1_ego[Csrs_borrow] == Links_borrow$threat1_alter[Csrs_borrow], 
                                   Links_borrow$threat1_ego[Tsrs_borrow] == Links_borrow$threat1_alter[Tsrs_borrow])$p.value #sig at .1, seek out similarit

#### In terms of shared perception of importance?
treateffect_borrow[1,11] <- mean(Links_borrow$important_surr_refug_ego[Csrs_borrow] == Links_borrow$important_surr_refug_alter[Csrs_borrow], na.rm=T)
treateffect_borrow[2,11] <- mean(Links_borrow$important_surr_refug_ego[Tsrs_borrow] == Links_borrow$important_surr_refug_alter[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,11] <- t.test(Links_borrow$important_surr_refug_ego[Csrs_borrow] == Links_borrow$important_surr_refug_alter[Csrs_borrow],
                                   Links_borrow$important_surr_refug_ego[Tsrs_borrow] == Links_borrow$important_surr_refug_alter[Tsrs_borrow]  )$p.value #sig at .5


#### In terms of shared positive view on topic?
treateffect_borrow[1,12] <- mean(Links_borrow$posview_homoph[Csrs_borrow], na.rm=T)
treateffect_borrow[2,12] <- mean(Links_borrow$posview_homoph[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,12] <- t.test(Links_borrow$posview_homoph[Csrs_borrow], Links_borrow$posview_homoph[Tsrs_borrow])$p.value

#### In terms of shared high interest in topic?
treateffect_borrow[1,13] <- mean(Links_borrow$highimp_homoph[Csrs_borrow], na.rm=T)
treateffect_borrow[2,13] <- mean(Links_borrow$highimp_homoph[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,13] <- t.test(Links_borrow$highimp_homoph[Csrs_borrow], Links_borrow$highimp_homoph[Tsrs_borrow])$p.value

#### In terms of mutliplexity?
treateffect_borrow[1,14] <- mean(Links_borrow$multiplicity_indicator[Csrs_borrow], na.rm=T)
treateffect_borrow[2,14] <- mean(Links_borrow$multiplicity_indicator[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,14] <- t.test(Links_borrow$multiplicity_indicator[Csrs_borrow], Links_borrow$multiplicity_indicator[Tsrs_borrow])$p.value


#### In terms of a count of layers?
treateffect_borrow[1,15] <- mean(Links_borrow$multiplicity_count[Csrs_borrow], na.rm=T)
treateffect_borrow[2,15] <- mean(Links_borrow$multiplicity_count[Tsrs_borrow], na.rm=T)
treateffect_borrow[3,15] <- t.test(Links_borrow$multiplicity_count[Csrs_borrow], Links_borrow$multiplicity_count[Tsrs_borrow])$p.value

colnames(treateffect_borrow) <- c("Links", "Alter was ref", "Alter knows ref", "Alter farmer", "Alter's views", "Alter's interest", "Relig homoph", "Language homoph", 
                                 "Refugee status homoph", "Refugee views homoph", "Interest homoph", "Positive views homoph", "High interest homoph",
                                 "In Multiple Layers", "Count of Layers")

rownames(treateffect_borrow) <- c("Control", "Treatment", "p-value")



## And build Appendix Table A6
xtable(t(treateffect_borrow), digits = 2)





#######
## And a histogram of the number of discussion partners
## Appendix Figure A1

plot(hist(D$numspokeref))



########
## Explore the overlap in the layers

## 1) Show the count of links in each layer that also appear in the others,
## maybe as a heatmap

## 2) Generate rows for table 6 that compare links used to those not used
## and verify with QAP


names(Links)
dim(Links)
dim(L)
table(Links$in_union)
table(L$in_union)

## First, let's create a multiplexity measure for our links
## L is our dataframe of links that are in the union social network

L$multiplicity_count <- L$in_meal + L$in_visit + L$in_rumor + L$in_borrow
L$multiplicity_indicator <- ifelse(L$multiplicity_count > 1, 1, 0)

table(L$multiplicity_count)

L %>% ggplot(aes(x = multiplicity_count)) +
  geom_histogram(color = "black" , fill = "grey65", binwidth =1) +
  theme_minimal()+
  labs(title = "Multiplexity", 
       x = "Number of layers in which a link appears", 
       y = "Count of Links") 

## Next, let's try building a stacked barplot showing the composition of the
## three layers


L <- L %>%
  mutate(multiplicity_comp = case_when(
    multiplicity_count == 1 & in_meal == 1 ~ "M",
    multiplicity_count == 1 & in_visit == 1 ~ "V",
    multiplicity_count == 1 & in_rumor == 1 ~ "R",
    multiplicity_count == 1 & in_borrow == 1 ~ "B",
    multiplicity_count == 2 & in_meal == 1 & in_visit == 1 ~ "MV",
    multiplicity_count == 2 & in_meal == 1 & in_rumor == 1 ~ "MR",
    multiplicity_count == 2 & in_meal == 1 & in_borrow == 1 ~ "MB",
    multiplicity_count == 2 & in_visit == 1 & in_rumor == 1 ~ "VR",
    multiplicity_count == 2 & in_visit == 1 & in_borrow == 1 ~ "VB",
    multiplicity_count == 2 & in_borrow == 1 & in_rumor == 1 ~ "BR",
    multiplicity_count == 3 & in_meal == 1 & in_visit == 1 & in_borrow == 1 ~ "MVB",
    multiplicity_count == 3 & in_meal == 1 & in_visit == 1 & in_rumor == 1 ~ "MVR",
    multiplicity_count == 3 & in_meal == 1 & in_borrow == 1 & in_rumor == 1 ~ "MBR",
    multiplicity_count == 3 & in_visit == 1 & in_borrow == 1 & in_rumor == 1 ~ "VBR",
    multiplicity_count == 4 ~ "MVBR"))

L %>%
  ggplot(aes(x = multiplicity_count, fill = multiplicity_comp)) +
  geom_histogram(binwidth = 1) +
  geom_text(aes(label = multiplicity_comp))

L_hist <- L %>%
  group_by(multiplicity_count, multiplicity_comp) %>%
  summarize(n = n())

L_hist %>%
  ggplot(aes(x = multiplicity_count, y = n, fill = multiplicity_comp)) +
  geom_col(color = "white") +
  geom_text(aes(label = multiplicity_comp), size = 3, position = position_stack(vjust = 0.5))+
  theme_minimal() + 
  labs(title = "Multiplexity", 
       x = "Number of layers in which a link appears", 
       y = "Count of Links") +
  guides(fill = "none")

## save as multiplexity.pdf

## Now, can we show the links used versus not separated by these?


L_hist2 <- L %>%
  mutate(in_spokeref = case_when(
    in_spokeref == 0 ~ "Not Used",
    in_spokeref == 1 ~ "Used")
  ) %>%
  group_by(in_spokeref, multiplicity_comp) %>%
  summarize(n = n())

L_hist2 %>%
  ggplot(aes(x = in_spokeref, y = n, fill = multiplicity_comp)) +
  geom_col(color = "white") +
  #geom_text(aes(label = multiplicity_comp), size = 3, position = position_stack(vjust = 0.5))+
  theme_minimal() + 
  labs(title = "Multiplexity", 
       x = "", 
       y = "Count of Links")+
  guides(fill = guide_legend(title = "Layers")) 
  
#save as multiplexity_byuse.pdf


summary(lm(in_spokeref ~ multiplicity_comp, data = L))

## And having the omitted category be any one:

L$multiplicity_comp_2plus <- L$multiplicity_comp
L$multiplicity_comp_2plus[which(L$multiplicity_count == 1)] <- "One"

L$multiplicity_comp_2plus <- as.factor(L$multiplicity_comp_2plus)
L$multiplicity_comp_2plus <- relevel(L$multiplicity_comp_2plus, ref = "One")







mult1 <- lm(in_spokeref~multiplicity_comp_2plus, data = L)

## Make a coefficient plot that orders the coefficients by their size
tidy(mult1, conf.int = TRUE) %>% 
  filter(term != "(Intercept)") %>%
  mutate(term = fct_reorder(term, estimate)) %>%
  ggplot(mapping = aes(x = estimate, y = term, xmin = conf.low, xmax = conf.high)) +
  geom_pointrange() +
  labs(title = "Layer Composition Predicting Use of Link to Discuss Refugees",
       x = "Coefficient Estimate from OLS, Omitted Category is Any Single Layer",
       y = "Layers") +
#  scale_y_discrete(labels = labs) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = .8)+
  theme_minimal()

## And now make a vector of variable labels, remembering that they 
## display from the bottom up

labs <- c("MB", "MVR", "MR", "BR", "MV", "VB", "VR", "MBR", "MVB", "MVBR", "VBR")

tidy(mult1, conf.int = TRUE) %>% 
  filter(term != "(Intercept)") %>%
  mutate(term = fct_reorder(term, estimate)) %>%
  ggplot(mapping = aes(x = estimate, y = term, xmin = conf.low, xmax = conf.high)) +
  geom_pointrange() +
  labs(title = "Layer Composition Predicting Use of Link to Discuss Refugees",
       x = "Coefficient Estimate from OLS, Omitted Category is Single-Layer Links",
       y = "Layers") +
  scale_y_discrete(labels = labs) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = .8)+
  theme_minimal()

## Save as multiplicity_predicting.pdf

tidy(mult1, conf.int = TRUE) %>% 
  filter(term != "(Intercept)") %>%
  mutate(term = fct_reorder(term, estimate)) %>%
  ggplot(mapping = aes(x = estimate, y = term, xmin = conf.low, xmax = conf.high)) +
  geom_point() + #this one turns off SEs
  labs(title = "Layer Composition Predicting Use of Link to Discuss Refugees",
       x = "Coefficient Estimate from OLS, Omitted Category is Single-Layer Links",
       y = "Layers") +
  scale_y_discrete(labels = labs) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = .8)+
  theme_minimal()
## Save as multiplicity_predicting_no_ses.pdf

## Make a table of the proportion of links in one layer that also appear in another

overlap_mat <- matrix(data = NA, nrow = 4, ncol = 4)
rownames(overlap_mat) <- c("Meal", "Visit", "Rumor", "Borrow")
colnames(overlap_mat) <- c("Meal", "Visit", "Rumor", "Borrow")

overlap_mat[1,2] <- sum((L$in_meal & L$in_visit) == 1) / sum(L$in_meal == 1) #prop meal in visit
overlap_mat[1,3] <- sum((L$in_meal & L$in_rumor) == 1) / sum(L$in_meal == 1) #prop meal in rumor
overlap_mat[1,4] <- sum((L$in_meal & L$in_borrow) == 1) / sum(L$in_meal == 1) #prop meal in borrow

overlap_mat[2,1] <- sum((L$in_visit & L$in_meal) == 1) / sum(L$in_visit == 1) #prop visit in meal
overlap_mat[2,3] <- sum((L$in_visit & L$in_rumor) == 1) / sum(L$in_visit == 1) #prop visit in rumor
overlap_mat[2,4] <- sum((L$in_visit & L$in_borrow) == 1) / sum(L$in_visit == 1) #prop visit in borrow

overlap_mat[3,1] <- sum((L$in_rumor & L$in_meal) == 1) / sum(L$in_rumor == 1) #prop rumor in meal
overlap_mat[3,2] <- sum((L$in_rumor & L$in_visit) == 1) / sum(L$in_rumor == 1) #prop rumor in visit
overlap_mat[3,4] <- sum((L$in_rumor & L$in_borrow) == 1) / sum(L$in_rumor == 1) #prop rumor in borrow

overlap_mat[4,1] <- sum((L$in_borrow & L$in_meal) == 1)/ sum(L$in_borrow == 1) # prop borrow in meal
overlap_mat[4,2] <- sum((L$in_borrow & L$in_visit) == 1)/ sum(L$in_borrow == 1) # prop borrow in visit
overlap_mat[4,3] <- sum((L$in_borrow & L$in_rumor) == 1)/ sum(L$in_borrow == 1) # prop borrow in rumor

## Let's can add in the 1s on the diagonal

overlap_mat[1,1] <- 1
overlap_mat[2,2] <- 1
overlap_mat[3,3] <- 1
overlap_mat[4,4] <- 1

## And output this as a table

xtable(overlap_mat, digits = 2)

#####
## And now let's assess the case for multiplicity mattering for link use

summary(lm(in_spokeref ~ multiplicity_count, data = L))
summary(lm(in_spokeref ~ multiplicity_indicator, data = L))

summary(lm(multiplicity_indicator ~ in_spokeref, data = L))

## Let's perform QAP to get these standard errors right


nsims <- 1000
set.seed(20240421)

numvlgs <- length(unique(D$village_number))

sim_reg_coefs_4 <- matrix(data = NA, nrow = nsims, ncol = 2)
colnames(sim_reg_coefs_4) <- c("multiplicity_count", "multiplicity_indicator")

L_sm2<- L %>% select(ego_hh, alter_hh, village_number, in_spokeref, multiplicity_count,
                    multiplicity_indicator)

## These are features of a node's links, not just of a node.  So now we need to shuffle
## these link attributes among the links, not among the nodes.


for (i in 1:nsims){
  
  ## within village, randomize the link attribute multiplicity count, shufflfed over the links
  
  L_sim <- L_sm2
  
  for(j in 1:numvlgs){
    
    ## Identify the rows of the attribute matrix that are part of village j
    thevlg <- which(L_sim$village_number == j)
    
    ## Shuffle these rows by using an index randomly shuffled
    
    shuffled_vlg <- sample(1:length(thevlg), length(thevlg), replace = FALSE)
    
    ## Now replace each of that village's link attributes with the same that have been shuffled
    
    ## Multplicity_count:
    L_sim$multiplicity_count[thevlg] <- L_sim$multiplicity_count[thevlg][shuffled_vlg]

  
  }
  
  ## Now L_sim contains a shuffled multiplicity_count
  ## Let's build the multiplicity indicator too
  
  L_sim$multiplicity_indicator <- ifelse(L_sim$multiplicity_count > 1, 1, 0)
  
  
  ## Now run the two simple regressions and store the coefficients from each
  sim_reg_coefs_4[i,1] <- coef(lm(multiplicity_count ~ in_spokeref, data = L_sim))[2]
  sim_reg_coefs_4[i,2] <- coef(lm(multiplicity_indicator ~ in_spokeref, data = L_sim))[2]

  
  if(i %% 100 == 0) print(i)
  
}

## Next, store the correct reference coefficient



coef_mult_count <- coef(lm(multiplicity_count ~ in_spokeref, data = L))[2]
coef_mult_indicator <- coef(lm(multiplicity_indicator ~ in_spokeref, data = L))[2]



## And now assess significance
sum(sim_reg_coefs_4[,1] >= coef_mult_count) #0
sum(sim_reg_coefs_4[,2] >= coef_mult_indicator)  #0

## And we can generate the two-sided QAP scores with:

sum(abs(sim_reg_coefs_4[,1]) >= abs(coef_mult_count))/1000 # 0
sum(abs(sim_reg_coefs_4[,2]) >= abs(coef_mult_indicator))/1000 # 0
head(sim_reg_coefs_4)

## And plot these

mult_regs <- as.data.frame(sim_reg_coefs_4)

mult_regs %>%
  ggplot()+
  geom_density(aes(multiplicity_count), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_mult_count, linetype = "dashed")+
  labs(title = "QAP, Multiplicity Count", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_mult1.pdf

mult_regs %>%
  ggplot()+
  geom_density(aes(multiplicity_indicator), alpha = .8, fill = "grey85")+ ## can add fill argument if we want a color
  geom_vline(xintercept = coef_mult_indicator, linetype = "dashed")+
  labs(title = "QAP, Multiplicity Indicator", 
       x = "Difference in Means from Permuted Data", 
       y = "") +
  theme_minimal()
#save as qap_reg2.pdf

## And build a table in the paper (Table 3):

mult_mat <- matrix(data = NA, nrow = 2, ncol = 4)
mult_mat[1,1] <- mean(L$multiplicity_count[srsL], data = NA)
mult_mat[1,2] <- mean(L$multiplicity_count[notsrsL], data = NA)
mult_mat[1,3] <- t.test(L$multiplicity_count[srsL], L$multiplicity_count[notsrsL])$p.value
mult_mat[1,4] <- 0 #calculated above


mult_mat[2,1] <- mean(L$multiplicity_indicator[srsL], data = NA)
mult_mat[2,2] <- mean(L$multiplicity_indicator[notsrsL], data = NA)
mult_mat[2,3] <- t.test(L$multiplicity_indicator[srsL], L$multiplicity_indicator[notsrsL])$p.value
mult_mat[2,4] <- 0 #calculated above

mult_mat <- as.data.frame(mult_mat)

xtable(mult_mat, digits = 2)
























